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ABSTRACT 

We present local two-dimensional (2D) and three-dimensional (3D) hybrid numerical simulations 
of particles and gas in the midplane of protoplanetary disks (PPDs) using the Athena code. The 
particles are coupled to gas aerodynamically, with particle-to-gas feedback included. Magnetorota- 
tional turbulence is ignored as an approximation for the dead zone of PPDs, and we ignore particle 
self-gravity to study the precursor of planctesimal formation. Our simulations include a wide size 
distribution of particles, ranging from strongly coupled particles with dimensionless stopping time 
t s = 51t st0 p = 1CP 4 to marginally coupled ones with t s = 1 (where Q is the orbital frequency, i s top is 
the particle friction time), and a wide range of solid abundances. Our main results arc: 1. Particles 
with t s > 10~ 2 actively participate in the streaming instability, generate turbulence and maintain the 
height of the particle layer before Kelvin-Helmholtz instability is triggered. 2. Strong particle clump- 
ing as a consequence of the streaming instability occurs when a substantial fraction of the solids arc 
large (r s > 10~ 2 ) and when height-integrated solid to gas mass ratio Z is super-solar. We construct a 
toy model to offer an explanation. 3. The radial drift velocity is reduced relative to the conventional 
Nakagawa-Sckiya-Hayashi (NSH) model, especially at high Z. Small particles may drift outward. We 
derive a generalized NSH equilibrium solution for multiple particle species which fits our results very 
well. 4. Collision velocity between particles with t s > 10~ 2 is dominated by differential radial drift, 
and is strongly reduced at larger Z. This is also captured by the multi-species NSH solution. Various 
implications for planetesimal formation are discussed. In particular, we show there exist two positive 
feedback loops with respect to the enrichment of local disk solid abundance and grain growth. All 
these effects promote planetesimal formation. 

Subject headings: diffusion — hydrodynamics — instabilities — planetary systems: protoplanetary 
disks — planets and satellites: formation — turbulence 



1. INTRODUCTION 

Planets are believed to be formed out of dust grains 
that collide and accrete into larger and larg er bodies 
in the gaseous protoplanetar y disks (PPDs) (|Safronovl 
119691: IChiang fe Youdinl 12009ft . The remarkable growth 
of dust into planets covers 40 orders of magnitude in 
mass, and can be divided into three regimes. At cen- 
timeter size or less, chemical bond and electrostatic 
forces allow small dust grai ns to stick to each other 
to form larger aggregate s (jDominik fe Tieleiisl Il997t 
iBlum fe Wuridl2000ll2008j) . At kilometer or larger sizes 
(i.e., planetesimals and larger bodies), gravity is strong 
enough to retain collision fragment s, leading to the for- 
mation of planetary embry o s /cor e s (I Wetherill fe Stewart! 
1981 iLissauer fe Stewart! fl993t iKokubo fc Idal 119981: 
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20041). and u l timat e ly to terrestrial an d 



giant planets dPollack et al.1 Il996t llda fc Linl [2004al lbl: 
iKenvon fc Bromlevl 120061 ) . The intermediate size range 
lies in the regime of planetesimal formation. This is prob- 
ably the least understood process in planet formation, 
largely because of solid growth in this regime is subject 
to a bottleneck known as the "meter size barrier" . 

In the intermediate size range, aerodynamic coupling 
between gas and solids is important. The gaseous disk 
is partially supported by a radial pressure gradient, and 



rotates at sub-Keplerian velocity, while solid bodies tend 
to orbit at Keplerian velocity. Consequently, solid bod- 
ies feel a headwind and drift radially inwards due to gas 
drag. The infall ti me scale is of the order 10 2 years for 
meter-sized bodies (|Weidenschillingl ri977T ). which poses 
strong constraint on the timescale of planetesimal for- 
mation. Moreover, the collision velocity between meter 
sized boulders and other bodies is large enough to re- 
sult in bouncing^ or fragmentation (jGiittler et al.l [20091 : 
IZsom et al.ll2010h . rather than growth. To overcome the 
meter size barrier, collective effects that form planetesi- 
mals out of meter size d or smaller bodies appear to be 
essential. For example, ICuzzi et al.l (|2001[ ) proposed the 
turbulent concentration of chondrule sized particulates 
by factors of up to 10 5 by extrapolating experimental re- 
sults to high Reynolds numbers. It such dense regions, 
mutual gravity of the particulates as a whole can over- 
come ram pressure and draw them together to form plan- 
etesimals (|Cuzzi et al.ll2008l ) , although the intermittency 
in th e turbulence might w ork against particle concentra- 
tion (lYoudin fc Shul^OOl . 

One favorable model of planetesimal formation in- 
volves gravitational instability (GI) in the settled 
dust layer in the midp lane of PPDs (ISafronovl 119691 : 
iGoldreich fc W ard 1973J). In the absence of turbulence in 
the disk, the dust layer would become thinner and thin- 
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mals by gravitational collap se and fragmentation. How- 
ever, as first pointed out by IWeidenschillingl (|1980t ). tur- 
bulence generated by vertical shear across the midplane 
dust layer (via the Kelvin-Hclmholtz instability, here- 
after KHI) prevents dust grains from continuously set- 
tling well before GI is able to operate. Based on the clas- 
sical criterion for the onset of the KHI and solar mctalic- 
ity for height-integrated dust to gas mass ratio (hereafter, 
solid abundance, denoted by Z), the maximum solid den- 
sity in disk midplane was found to be generally 1-2 orders 
of ma gnitude lower than the Roche dens ity for the onset 
of GI IjSekivalllQM lYoudin fc Shul EfflOlh 1 . Inclusion of 
Coriol is force ([Gomez fe Ostrikerl|2005f) as well as radial 
shear ([Chi ang 2008; Ba rrancoll2009|) do not alter the con- 
clusion qualitatively. It appears that increasing the local 
solid abundance by a factor of 2-10 times solar is needed 
for this mechanism to operate 2 . Thi s factor may be 
achievable by photoevaporation of gas ([Throop fc Ballvl 
120051: lAlexander fc Armitagell2007fl . and by the radial 
variations of orbital drift speeds induced b y gas drag 
([Youdin fc Shul[200l lYoudin fc Chiang|l2004fl . 

An important ingredient of particle-gas interaction in 
the midplane solid layer is the backreaction from par- 
ticles to the gas. The momentum feedback from solids 
to gas is responsible for KHI which tends to maintain a 
finite thickness of the solid layer. When the solids are 
not too strongly coupled to the g as, the backreaction 
leads to a powerful drag instability ([Goodman fc Pindorl 
l2000f ) . n ow termed the "streamin g instability" (here- 
after SI. lYoudin fc Goodman! [20051 ). The most remark- 
able feature of the SI is that it ver y efficiently concen- 
trates particles into dense clum ps ([Youdin fc Johansenl 
l2007t [johansen fc Youdinll2007t ). and enhances local par- 
ticle density by a factor of up to 10 3 . Such enhance- 
men t in particle density is su fficient to trigger GI, 
and iJohansen et al.l (J20071 . l2009f ) found in their simula- 
tions that planetesimals form rapidly once self-gravity is 
turned on. The sizes of the planetesimals formed in the 
simulations are about a few hundreds kilometers, consis- 
tent with constraints deduced from observations of as- 
teroid and Kuiper belt objects th at planetesimals are 
formed big ([Morbidelli et alj|2009h . These results pro- 
vide a very promising path for forming planetesimals by 
SI followed by gravitational collapse. 

Planetesimal formation is also affected by external 
turbulence in PPDs. The typical ma ss accretion rate 
of_10" 8±1 MQ yr _1 for T-Tauri stars ([Hartmann et al.l 
ll998[ ) indicates efficient angular momentum transport in 
PPDs. Magnetic field seems certain to play a crucial 
role in the transport process, most noticeably by the 
magnetorotational instability (MRI) ([Balbus fc Hawleyl 
119911: lHawlev fc Balbuslll99lT ). The turbulence gener- 
ated by MRI strongly affect the s ettling of sma l l dust 
grains ([Fromang fc Nelson! 120091: iBalsara et all 120091: 

1 The Roche density criterion for the onset of GI may not apply 
to the dust sublayer due to the drag interaction between gas and 
solids, and |Youdin (2005a b]) showed that GI can occur at lower 
densities with smaller growth rate, although turbulent diffusion of 
solids is ignored in his calculation. 

2 See also the most recent results by[Lec ct al. (2010b]) who stud- 
ied the onset of KHI from more realistic dust density profiles from 
dust settling. 



iTillev et al.l l2009t) . but more interestingly, it promotes 
the concentration of dec i meter to meter sized bodie s 
([Fromang fc Nelson! [20051: IJohansen et al.ll2006bl l2007t) . 
PPDs arc, however, only weakly ionized. The main ion- 
ization sources such as cosmic rays and X-rays from 
the protostar only ionize the surface of the disk, mak- 
ing the surface layers "active" to MRI driven turbu- 
lence, w hile the midpla ne remains poorly ionized and 
"dead" (jGammid Il996[ ). Accretion is therefore lay- 
ered and mainly proceeds in the active zone. More- 
over, the presence of small dust grains substantially in- 
creases disk r esistivity and reduces the extent of the 
active layer (ISano et all 120001: lllgner fc Nelsonl 120061 : 
iSalmeron fc Wardldl2008t iBai fc Goodmanll2009ri . These 
non-ideal MHD effects due to partial ionization and dust 
resistivity, as well as the layered accretion structure in 
PPDs tremendously complicate the story of planetesimal 
formation. 

In this paper, we consider a local patch of PPDs 
and study the dynamics of gas and solids in the disk 
midplane. We perform shearing box hybrid simula- 
tions with both ga s and particles using the Athena code 
([Stone et al.l 120081 ) . The implementation of the parti- 
cle mod ule and code tests are presented in lBai fc Stone! 
(|2010al) . The inclusion of backreaction from particles to 
gas allows us to investigate both the SI and KHI simul- 
taneously. The local model is necessary for studying SI 
because the scale of particle clumping is much smaller 
than gas scale hei ght and requires at l east 16 cells to be 
properly resolved ([Bai fc Stone!l2010al ). The self-gravity 
from particles is neglected. Although self-gravity will 
ultimately play an important role in planetesimal forma- 
tion, our focus is its precursor: clumping of particles. 
Neglecting self-gravity also has the advantage that our 
results can be easily scaled to different disk parameters 
and have very broad applications (see §2.2p . We have 
also neglected the thermodynamics in our work, which 
may affect the buoyancy of the gas, but the dynamics 
of th e particles are generally unaffected ([Garaud fc Lin! 
l200l . 

Our ultimate goal is to build the most realistic local 
model of PPDs possible, including all of the non-ideal 
MHD effects as well as dust grains/solid bodies in a self- 
consistent manner. In this paper, however, we focus on 
the dynamics in the dead zone, and therefore can neglect 
MHD. This simplification is justified in two ways. First, 
conductivity calculations have shown that the inner part 
of PPDs (r < 10AU) almost always contains a dead zone 
(|Bai fc Goodman! 120091 ; [Turner fc Drake! l2009f >. Second, 
this approach separates the hydrodynamic effects (SI) 
from non-ideal magnctohydrodynamic (MHD) effects, 
which sets the foundation for more sophisticated work. 
In reality, the dynamics in the dead zo ne can be affected 
by th e turbulence in the active layer ([Fleming fc Stone! 
120031 ). For example, the gas motion in the disk mid- 
plane may exhibit strong low-frequency (compared with 
orbital frequency fl) vertical oscillations excited by the 
turbulence in the upper layer, and no coher ent anti-cyclic 
vortices are found (|Oishi fc Mac Lowll2009D . Its influence 
to the dynamics of the solids is not clear and is left for 
future investigations. 

An important ingredient of our simulations is the size 
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distribution of particles. A wide size distribution of dust 
grains from micron to millimeter or centimeter size in 
the PPDs is routinely deduced from the modeling of 
their spectral energy d i stribution (SEP ) (IChiang et al.l 
120011 : iTesti et al.l 120031 : iD'Alessio etHI I2006D . Theo- 
rctical modeling of dust coagu lation also result in a 
broad range of particle sizes (|Dullemond fe Dominikl 
120051: iBrauer et al.|[2008at iBirnstiel et alll2010D . In the 
most recent work that incorporate s up-to-date labora- 
tory collision ex periment results (|Giittler et al.l 120091 : 
IZsom et al.ll2010l ). the particle size range that dominates 
the total solid mass spans about 1-3 magnitude, typi- 
cally from sub -millimeter to decimeter ran ge. We note 
that although iJohansen et al.l (|2007l l2009f) also consid- 
ered a size distribution of particles, their particle size is 
relatively large and the size range is narrow (maximum 
particle size is 4 times the smallest). In this paper, we 
choose the particle size range to span 1-3 orders of magni- 
tude, and we assume uniform particle mass distribution 
in logarithmic size bins. Our choice of the particle size 
distribution roughly agrees with outcome of coagulation 
model calculations and serves as a first approximation of 
reality. We perform a parameter survey on particle size 
range and height-integrated particle to gas mass ratio (or 
solid abundance) that cover a substantial fraction of pa- 
rameter space relevant to planetesimal formation. These 
simulations sclf-consistcntly include the mutual interac- 
tions between gas and parti cles of all sizes (e xtending 
the early analytical work by ICuzzi et al.1 119931 who as- 
sumed all particles are passive), and will help us better 
understand the environment and precursor of planetesi- 
mal formation. 

We perform both two-dimensional (2D) and three- 
dimensional (3D) simulations, where the 2D simu- 
lations are axisymmetric (i.e., in the radial- vertical 
plane). We note that K HI is most prominen t in 
the azimuthal- vertical plane (jJohansen et al.ll2006al ). al- 
though fully capturing K HI requires fully 3D simula- 
tions including ra dial shear (Chiang 2008; Bar rancol2009l : 
iLee et all l2010al) . On the other hand, 2D simulation 
in the radial-vertical plane is sufficient to cap t ure S I 
(|Youdin fc Goodman! 120051: IJohansen fc Youdinl I2007D . 
While 3D simulations are necessary to capture all possi- 
ble physical effects in the disk midplane layer, we show 
in £13.11 that KHI is unlikely to be present in all our 3D 
simulations, because the turbulence generated by SI is 
strong enough to prevent the particles from further set- 
tling to trigger KHI 3 . Therefore, 2D simulations are also 
a valid approach to the problem, and are much less time- 
consuming than the corresponding 3D runs. Moreover, 
comparison between 2D and 3D simulations can be used 
for discerning multi-dimension effects, and as a guidance 
for future studies. 

This paper is organized as follows. In Sj2j we describe 
our simulation method, model parameters and scaling 
relations. We also describe the basic properties of the 
saturated state in all our simulations. We study var- 
ious aspects of our simulations in the subsequent four 
sections. In SJ3] we discuss the vertical structure of the 

3 This is no longer true if all particles are strongly coupled to 
gas, in which case the SI is much weaker. 



particle layer. In particular, we address the question of 
what is the dominant process of the midplane dynamics, 
KHI or SI? We further analyze which particles are ac- 
tively participating in the instabilities, and which parti- 
cles behave only passively. In Sj4] we study the conditions 
for forming dense clumps from the SI, which preludes 
planetesimal formation. The composition and dynamics 
of the dense clumps is also analyzed. $5] deals with the 
radial transport of particles, including both radial drift 
and radial diffusion. We study particle collision velocities 
in ij6] We conclude our paper in Sj7] by summarizing our 
results and discussing various implications for planetes- 
imal formation. In particular, we summarize the logical 
connections between various physical effects that may en- 
hance each other and promote planetesimal formation. 

2. METHOD AND SIMULATIONS 

2.1. Formalism 

We consider local PPD models and formulate the equa- 
tions of gas and solids using the shearin g sheet approx- 
imation (|Goldreich fc Lvnden-Beillll965f ). We choose a 
local reference frame located at a fiducial radius, coro- 
tating at the Kcplcrian angular velocity 57. The dynam- 
ical equations are written using Cartesian coordinates, 
with x,y, z denoting unit vectors pointing to the radial, 
azimuthal and vertical direction, where $1 is along the 
z direction. The gas density and velocities are denoted 
by p g ,u in this non-inertial frame. We include a dis- 
tribution of particles coupled with gas via aerodynamic 
drag, where the velocity of particle i is denoted by Vi. 
The drag force is characterized by stopping time £ s top, 
and equals (u — v)/t s top per unit particle mass. Particles 
with different sizes have different stopping times, labeled 
by subscript "fc" . Back reaction from the particles to gas 
is included, which is necessary for the study of KHI and 
SI. In this non-inertial frame, the equations for the gas 
read 



^+V-( Pg u) = 



^ + V.( PgUU + PI) 



Pg 



2uxCl + il 2 xx - Vt 2 zz + y^ 



(k 



v k -u 



stop.fc 



(1) 



(2) 



where the source terms include Coriolis force, radial tidal 
potential as well as disk vertical gravity. The last term 
in the momentum equation represents the backreaction 
(or momentum feedback) from particles to gas: tk and 
v~k denote the local mass density and velocity of particles 
of type fc. In this paper we neglect the effect of magnetic 
fields and focus on the inter action between gas and solids 
in the dead zone of PPDs (iGam mie 199(f). An isother- 
mal equation of state for the gas is used throughout this 
paper, where P = p g c 2 s and c s is the isothermal sound 
speed. 

Similarly, the equation of motion for particle i of type 
fc can be written as 

dvi „ ,_ 2 „ -.q „ Vi — u 

—— = —2i]VkHx + 2vi x Si + HxiX — ilziZ . 

Ctfc tstop,A: 

(3) 
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In the above equation, we have added an inward force 
term —2rjVK^lx to mimic th e effect of an outwar d radial 
pressure gradient in the gas (|Bai fe Stonell2010al) . where 
tjvk is the difference between gas velocity and the Kep- 
lerian velocity in the absence of particles. This term will 
shift both gas and particle azimuthal velocities by tjvk 
relative to those in the real system. To avoid confusion, 
we always use u and v to denote velocities that corre- 
sponds to the real system (i.e., subtracting the azimuthal 
velocity component from the simulation by tjvk). Parti- 
cle self-gravity is ignored as we focus on the dynamics in 
the midplane of the PPD dead zone and precursor of the 
planetesimal formation. 

In our simulations, we have applied an or- 
bital advection algorithm for both gas and particles 
(jStone fc Gardiner! l20ltil iBai fc StoneTl2010aD . and the 
actual velocities used in the simulation are measured 
relative to the linearized Keplerian shear flow: u' = 
u + (3/2)ilxy for gas flow, and v\ = Vi + (3/2)rtxy for 
individual particles. 

2.2. Scaling Relations 

Measuring velocities in units of the sound speed, time 
in units of Q~ x , and length in units of the gas scale height 
H g = c s /fl, the parameters in the problem are reduced 
to the following: 

1. The dimensionlcss particle stopping time r^ = 
^stop.fc for particle species fc. 

2. The solid abundance parameter Z\. for each par- 
ticle species, which measures the height-integrated 
particle to gas mass ratio. 

3. The parameter characterizing the strength of the 
radial pressure gradient II = r\r j H g = tjvk /c s . 

Below, we apply a disk model and provide the scaling 
relation between the disk model parameters and these 
dimensionless parameters used in our simulation. 

We adopt a generalized solar nebular model where 
the disk is vertically isothermal and all the disk quan- 
tities have a powe r law dependence on the radius 
(|Youdinfc Shul 12001 



1700/ 9 rl£ g cm" 2 



(4) 



T = 280 f T rj£ K , 
M, = f M M Q . 
where E g is the gas surface mass density, T is the disk 
temperature, M* is the mass of the central star, and 
tau = r/lAU. These parameters fix the disk model. Al- 
though the global disk profile may not follow the simple 
power law form, we can always approximate a local patch 
of the disk in the above form, which is very general. 
In the standard mini mum-mass solar nebular (MMSN) 
model (|Havashilll98l . we have b = 3/2, c = 1/2, f T = 
f g = /m = 1. The radial profiles of other physical quan- 
tities are 



n = 2tt/ 



V2 -3/2 1 



M 'AU 



vr 



VK 



30/ Af r ATT km s 



'M 'AU 

- f ' r7°/ km s~ l 
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where in the calculation of the sound speed, we assume 
the mean molecular weight p = 2.33. 
The background gas density profile is 

p g , b (r, z) = -J^- exp(-z 2 /2H 2 g ) , (6) 



2irH g 

where subscript "6" denotes "background". Using this 
gas density and sound speed, one can derive the radial 
pressure gradient in the gaseous disk, thus obtain the 
amount of reduction t\vk in the gas rotation velocity. 
After some algebra, we can derive the pressure length 
scale parameter 

ldlnP c s (3 + 2b + c 3 - c z 2 \ c s 

2 dlnr v K ~ V 4 4 H 2 g J v K (J) 



0.054/^ 2 /m 1/2 ^ /4 



AU 



Note that n = II(r, z) depends on both radius and 
height. Nevertheless, in this paper, our simulation box is 
concentrated in the disk midplane where z <C H g , there- 
fore we can neglect the dependence of II on z. In the 
last equation of the above formula, we have applied the 
power law indices of the MMSN model. The dependence 
on disk temperature /t, stellar mass Jm as well as disk 
radius r is relatively weak. It is worth mentioning that 
the dependence of II on disk mass is only through the 
surface density profile parameter 6, free from the scal- 
ing parameter f g . Therefore, the value II rj 0.05 should 
apply to a wide range of disk models. 

Next we consider the scaling relations for the dimen- 
sionless stopping time. Because the gas motion in PPDs 
is expected to be subsonic, the relevant drag laws from 
the gas to the so lids in PPDs are the Epstein drag 
law (|Epsteinlll924D . which applies when particle size is 
smaller than the gas mean free path, and the Stokes drag 
law, which applies for larger bodies. We assume all solid 
bodies have spherical shapes, then the stopping time in 
these two regimes can be expressed as dWeidenschillind 
[1971 



L stop 



p s a 

Pg f S \ 
4:p s a 2 

9p g C s \r, 



a < 9A m /4 (Epstein regime), 
a > 9A m /4 (Stokes regime). 



(8) 

where p s s=a 3g cm 3 and a are the density and radius of 
the solid body, A m = {n g a)~ l = pmn/ p g o~ is the mean 
free path of the gas, and a ~ 2 x 10~ 15 cm 2 is the molec - 
ular collision cross section (jChapman fc Cowhnd[l970() . 
From the above equations, we see that the particle stop- 
ping time depends linearly on gas density in the Epstein 
regime. Nevertheless, the gas density can be regarded as 
constant near the disk midplane where we study. There- 
fore, in our local simulations, we can safely take i s top as 
depending on particle size a only. 

To better handle the relation between particle size and 
its corresponding stopping time, we express the relation 
between r s = f2t s t op and a by applying our disk model. 
The result is 



r, = max 



H g = 3.4 x 10-Vt'Vm ' 



-2 f l/2 f -l/2J3-c)/2 A[] _ 



4.4 x 10- 3 o cm /- 1 rl u 



1 A v in- 3 ^ 2 f _1 /2j; „(c-3)/2 

1.4 x 10 a cm / T jM r AU 



(9) 
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where a cm is the particle radius measure in centimeter. 
In the MMSN model, at 1 AU, particles smaller than 3cm 
are in the Epstein regime. At larger radii, the Epstein 
regime applies to much larger particles. 

2.3. Simulation Setup 

Fiducially, we consider the MMSN model at 1 AU, and 
set the pressure length scale parameter n = 0.05. This 
parameter is kept fixed in all our simulations. Instead of 
considering a particle size distribution in radius, we con- 
sider the distribution in r s . Then one can easily translate 
it into particle radius given the parameters of the disk 
model. We discretize a continuous particle size distri- 
bution into a number of bins. Each bin covers half a 
dex in r s in the logarithmic scale. For simplicity, we 
assume a uniform particle mass distribution across the 
bins, that is, all the particle bins (or particle species) 
have equal amount of mass. The parameters for the size 
distribution is therefore the minimum and maximum di- 
mensionless stopping time r m i n and r max (translated to 
a m i n and a max respectively). Physically, our assumption 
means that most of the mass of the solids resides in the 
size range between a m ; n and a max and roughly follows a 
flat distribution in logarithmic scale. To control the total 
particle mass, we use the total solid abundance parame- 
ter 



Z 



'type 

E 
fe=i 



Zh 



with Zi. 



Z/N< 



type , 



(10) 



where At ypo is the number of particle types (bins) . Cur- 
rently the best estim ate of the solar metallicity is about 
0.015 (|Lod dcrs 2003J). A substantial fraction of the metal 
elements may reside in dust grains and grow into larger 
bodies. In our simulations, we consider three abundance 
values Z = 0.01,0.02 and 0.03. This choice covers a 
relatively wide range of disk mctallicities. Moreover, be- 
cause our simulation focuses on a local patch in a PPD, 
the local abundance may not necessarily be equal to the 
averaged value in the PPD. 

As we explained in SJT] we perform simulations in both 
2D and 3D. Our 2D simulations are in the x-z plane 
(i.e. axisymmetric). Details of the implementation and 
co de tests of the parti cle-gas hybrid scheme are given 
in iBai fc Stone! (|2010a[) . Ou r simulations use th e stan- 
dard shearing box approach (|rlawlev et al.lH995l ). where 
the radial boundary condition is periodic with azimuthal 
shear. Azimuthal boundary conditions are periodic. Ver- 
tical gravity is included in our simulations, and we choose 
reflection boundary c ondition in the z dire ction, which 
is the same as that in lJohansen et all (|2009| ). In general, 
we use 256 cells in the ra dial (and azimuthal, if applica- 
ble) direction. Guided by IBai k, Stone! (|2010al ). properly 
resolving the SI with t s =0.1 requires about 128 cells 
per pressure length scale rjr. With this required resolu- 
tion, our simulation box size is typically small, spanning 
only about 2 — 4ryr. Such small box size is also necessary 
to capture the typical wavelength of the KHI, if present 
(|Johansen et al.ll2006aH . In our simulations, we generally 
use N p — 65536 particles per type for 2D simulations and 
N p = 3145728 particles per type in 3D runs (in which 
cases -/Vtypo = 7). Larger N p are used when N t y pc is 



smaller to keep the total number of particles similar in all 
our simulations. Our choice of particle number guaran- 
tees at least one particle per cell per particle type around 
the disk midplane, a s required for numerical convergence 
(|Bai fc Stonell2010a|K 

In our simulations, we set the initial particle den- 
sity profile to be a Gaussian centered on disk midplane 
with scale height Hd = 0.015iJ g for all particle types. 
The particle and gas velocities are computed from a 
multi-species Nakagawa-Sckiya-Hayashi (NSH) equilib- 
rium, where the classic al single-species NSH equilibrium 
(jNakagawa et al.lll986| ) solution is generalized to include 
multiple species of particles (see Appendix [A}. Note that 
different particles have different velocities, and the veloc- 
ities of particles and gas depend on z. 

The choice of our simulation box size and boundary 
conditions in the vertical direction merit further discus- 
sion. In the simulations, gas-particle interaction in the 
disk midplane generates turbulence and excites vertical 
motions in the gas. Ideally the vertical box size should 
extend to a few H a , similar to w hat is used for MRI 
simulations (e.g.. lStone et al.|[l996[ ). however, this would 
make 3D simulations too expensive. We have conducted 
a series of tests in 2D with a single particle species r s = 1 
using different vertical box sizes and either reflecting or 
periodic boundary conditions. In both cases, particles 
settle to the disk midplane with a spatial distribution 
reminiscent of sinusoidal waves that slowly drift in the 
radial direction. We find that the particle scale height is 
more intermittent when using periodic boundary condi- 
tions. Moreover, periodic boundary conditions appear to 
suppress asymmetric modes in the gas azimuthal veloc- 
ity around the disk midplane. Using reflecting boundary 
conditions, we find essentially no difference between the 
particle scale heights and clumping properties obtained 
from different vertical box sizes once the box height is 
much larger than the particle scale height, although it 
takes longer for the system to reach a quasi-steady state 
when a larger vertical box size is used. The drift veloc- 
ities of the wave-like pattern of particles do differ when 
different vertical box sizes are used, but they are unlikely 
to affect the properties discussed in $3] to $6] 4 . Guided 
by these results, as long as the vertical boundary of our 
simulation box is well above the scale height of all par- 
ticle species, one should get converged results from the 
simulations. 

Table [1] lists the parameters of all of our simulations. 
Our runs arc labeled using names with the form KxyZz- 
nD, where x, y are integers corresponding to r m - m = 
10 -a: , and T max = 10~ y , z = 100Z represents the solid 
abundance, and n = 2 (n = 3) denotes 2D (3D) simu- 
lations. When referring to simulations with fixed x and 
y but all possible values of z and/or n, we omit the Zz, 
and/or the nD, parts of the names. We focus on two 
groups of runs. In the first group, the maximum particle 
stopping time is r max = 0.1. We use 7 particle species to 
span three orders of magnitude in stopping time (down 
to r m j n = 10 -4 ) for the series of runs labeled R41, while 
in the series labeled R21, we use three particle species to 

4 Similar tests have been performed using the Pencil Code with 
the same conclusions (A. Johansen, private communication, 2009). 
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span one order of magnitude in stopping time (down to 
Tumi = 10~ 2 ). In the second group of runs, the maximum 
particle stopping time is r max = 1.0, and the minimum 
stopping time is chosen to be r m j n = 10~ 3 (R30) or 0.1 
(RIO). In each series of runs (R41, R21, R30, RIO), 
we perform three 2D simulations with Z = 0.01,0.02 
and 0.03, and two 3D simulations with Z = 0.01 and 
Z = 0.03. Because of a smaller r max in the first group, 
higher resolution is needed to resolve the SI. 



2.4. Simulation Runs and Saturation 

To determine when a saturated state is reached in each 
simulation, we monitor the particle vertical scale height 
H Pt k for each particle species fc, defined as the rms value 
of the z coordinate of all particles. Saturation occurs 
when particle settling and turbulent diffusion are in bal- 
ance, so that the scale height of all particle species is 
steady. In Figure Q] we show the time evolution of the 
vertical scale height for each particle species (marked by 
different colors) in all our runs. Solid and dashed curves 
represent 2D and 3D simulations respectively. We see 
that most of the 2D runs saturate within about 50 or- 
bits 5 . The 3D simulations are very time consuming, so 
we run them for shorter periods. From Figure [T] all 3D 
runs saturate before we terminate the simulations, al- 
though some just barely so. 

In the last column of Table [1] we provide the time of 
saturation T s (in parentheses) for each simulation. Un- 
less otherwise stated, we will perform data analysis in 
the time interval between the saturation time T s and 
the end time of the simulation T e . In the R41Z3-2D 
and R21Z3-2D runs, there are sudden jumps in par- 
ticle heights followed by settling, and this process re- 
peats over time quasi-periodically. Averaging over many 
cycles is required to reduce the influence of these in- 
termittent "bursts". The vertical distribution of the 
smallest particles in the R41-3D and R21-3D runs with 
Z = 0.03 are not fully saturated at the end of our sim- 
ulations. Nonetheless, the scale heights of the largest 
particles (which dominate the dynamics) in these runs 
have reached steady state, therefore we consider them to 
be saturated. 

Before presenting a detailed data analysis, we show the 
distribution of particles at the end of our simulations in 
Figure [2j Results from 3D runs are shown by projecting 
particle positions in three orthogonal directions. The 
number of particles plotted is much less than the actual 
number of particles used in the simulation. The trends 
in particle scale height evident in Figure Q] can be clearly 
seen: particles with small t s are diffused to larger heights. 
Note that we overplot larger particles on top of small 
particles, so that small particles near the midplane are 
less visible. The SI is present in all the simulations, and 
we will discuss various aspects of Figure[5]in the following 
sections. 

5 For run R41Z1-2D, the diffusion time of the smallest parti- 
cles with t s = 10 — 4 is very long and their H p still increases after 
1200C -1 . Nevertheless, the dynamics is dominated by the largest 
particles with r s > 10 -2 , and the scale heights of these particles 
has reached steady state. 



3. VERTICAL STRUCTURE OF THE DUSTY MIDPLANE 

LAYER 

3.1. Kelvin- Helmholtz Instability or Streaming 
Instability? 

The source of turbulence responsible for stirring up the 
particles can in principle be due to both KHI and SI. It is 
important to decipher which instability is the dominant 
process. Generally speaking, the onset of SI requires the 
averaged particle to gas mass ratio e > 1. The strength 
of the instability decreases as the averaged particle size 
becomes smaller, and vanishes asr s — > 0, for which the 
dust and gas behave as a single fluid. The onset of KHI 
requires a steep vertical profile of gas azimuthal velocity, 
which generally corresponds to larger dust to gas mass 
ratio at disk midplane. In our simulations, a substantial 
fraction of the particles have a relatively large stopping 
time with r s > 10~ 2 , and SI clearly plays an important 
role in the generation of disk midplane turbulence. It 
remains to study whether KHI is present and whether 
KHI is dynamically important. 

The classical result on the onset of KHI in a vertically 
stratified disk is based on the Richardson number criteria 
(|Chandrasekharlll961[ ) 



Ri 



x,y 



9_ {dp/dz) 
P (du x , y /dz) 



(11) 



where we define the Richardson number from radial and 
azimuthal velocity shear, as indicated by subscripts x,y. 
In the above equation, g = VI 2 z is the vertical gravita- 
tional acceleration, and p is the effective fluid density 
(see discussion below). The Richardson number mea- 
sures the amount of work required to overturn the fluid 
(numerator) in comparison to the amount of free energy 
available in the vertical shear (denominator). For Carte- 
sian flow with no rotation, the necessary condition for 
instability is given by Ri < Ri crlt = 1/4. This crite- 
ria no longer holds when rotation (Coriolis force) and 
radial shear (differential rotation) are included, espe- 
cially when the rotation frequency Q is comparable to the 
Brunt- Vaisala frequency of buoyant oscillations. Gen- 
erally speaking, the Cor iolis force distablizes the fluid 
(|G6mez fc Ostrikerll2005l ). while radial shear acts to sta- 
blize the fluid. iLee et al.l (|2010aD found that Ri CTit is 
typically smaller than 1/4 and is roughly proportional to 
dust to gas mass ratio at disk midplane. In this paper, we 
adopt the critic al Richa r dson number to be Ri crlt = 0.1 
as suggested by IChiand (|2008t ) . 

The Richardson number criterion is based on a single- 
fluid, in which case p simply represents fluid density. 
With the addition of perfectly coupled dust, the dust- 
gas system behaves as a single fluid, where the dust con- 
tributes to the mass but not the pressure of the fluid, 
thus p = p g + p p . When particles are not perfectly cou- 
pled, the definition of p becomes somewhat ambiguous, 
but we expect p g < p < p g + p v . Below we provide a sim- 
ple formula for p in this regime that reduces the above 
two limiting cases when p p = and when r s — > 0. 

In the absence of any turbulence and vertical grav- 
ity, the equilibrium state between gas and dust (with 
fixed stopping time) i s described by the NSH solution 
(|Nakagawa et al.lll986l ). In particular, the azimuthal gas 
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Run parameters. 








Run 


Z/0.01 1 


7m in 


^max 


JVtype 


J-Jx X -Lsy X -L/2 


N x x N y x N z 


N p 4 


T e (T s ) 5 


R41-2D 


1,2,3 


10~ 4 


io- 1 


7 


0.1 x - x 0.2 


256 x 1 x 512 


6.6 x IO 4 


1500(1200) 


R41-3D 


1,3 


10- 4 


io- 1 


7 


0.1 x 0.1 x 0.2 


256 x 256 x 512 


3.1 x IO 6 


250(200) 5 


R21-2D 


1,2,3 


IO" 2 


HT 1 


3 


0.1 x - x 0.2 


256 x 1 x 512 


9.8 x IO 4 


900(600) 


R21-3D 


1,3 


10- 2 


io- 1 


3 


0.1 x 0.1 x 0.2 


256 x 256 x 512 


6.3 x IO 6 


250(200) 5 


R30-2D 


1,2,3 


io- 3 


i 


7 


0.2 x - x 0.3 


256 x 1 x 384 


6.6 x IO 4 


1200(900) 


R30-3D 


1,3 


io- 3 


i 


7 


0.2 x 0.2 x 0.3 


256 x 256 x 384 


3.1 x IO 6 


450(300) 


R10-2D 


1,2,3 


io- 1 


i 


3 


0.2 x - x 0.3 


256 x 1 x 384 


9.8 x IO 4 


900(600) 


R10-3D 


1,3 


io- 1 


i 


3 


0.2 x 0.2 x 0.3 


256 x 256 x 384 


6.3 x IO 6 


450(300) 5 



runs with Z- 



1 Total particle to gas mass ratio, divided by 0.01. 

2 Number of particle species. 

3 Domain size, in unit of gas scale height H g = c s /Q. Note we have fixed II = r]r/H g = 0.05. 

4 Number of particles per species in the simulation box. 

5 Total run time in unit of f2 — 1 . The number in the parenthesis indicates the time of saturation. 
=0.03, we have T e = 280 and T a = 240. For run R10-3D, we have T e = 500 and T s = 450. 



For R41-3D and R21-3D 




200 400 600 800 200 400 600 800 200 400 600 800 200 400 600 800 
t(OT 1 ) t(£T 1 ) t(£T 1 ) t(£T 1 ) 



Fig. 1. — The time evolution of particle scale height H p for all simulations. Different colors represent different particle species, and 
particles with smaller t 3 have monotonically larger values of H p (see Table [T] for reference). Results from 2D simulations arc plotted with 
solid curves, while dashed curves show 3D results. Note that we run 2D simulations much longer than those in 3D, and the vertical scale 
in the top, middle and bottom panels are different. 
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Fig. 2. — The distribution of particles at the end of all our simulations. For each 3D run (shown in the leftmost and rightmost panels), 
we show the projected positions of a subset of particles in three orthogonal directions, while each 2D run is shown in one panel in the 
center. Different particle species are marked with different colors, and the color coding is the same as that used in Figure [l] Large red 
dots in a few plots (corresponding to the simulation runs that exhibit strong particle clumping) indicate the densest point in the particle 
clump. The unit of length in all panels is r\r. Note that the vertical size of our simulation box is larger than shown in this figure. 
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velocity relative to Keplerian velocity is given by 



1- 



(l + e) 2 +- 



rjv K ■ 



(12) 



where e = p v j p g , and prime means Keplerian velocity 
is subtracted. For convenience, we define Au y = —uL. 
For perfectly coupled particles, r s = and we find p g + 
Pp = PgV v K /Au y . For particles with finite stopping time, 
Au y becomes closer to tjvk, which reflects the fact that 
the particle-gas coupling is weaker so that gas velocity 
shifts towards the dust-free value. Therefore, Au y can 
be regarded as an indicator of particle-gas coupling. In 
this spirit, we define the effective gas density as 



„off _ 



Pg 



T]VK 

Au„ 



(13) 



It is trivial to check that in the limit e — >• 0, p — > p g . 
and p cfi — >• p g + p p when t s — >• 0. In the calculation of 
the Richardson number, we substitute p by p . Since p g 
is nearly constant over the height of our simulation box, 
equation (fTT|) becomes 



Ri 



x,y 



fl 2 z (duy/dz) 
Au y (du x .yjdzf 



(14) 



where the overbar means averaging over the horizontal 
plane. Note that Ri depends on z. 

Before calculating the Richardson number profile from 
our simulations, we first return to the spatial distribu- 
tion of particles in Figure [2] In 2D simulations, we 
see that the distribution of particles around the disk 
midplanc is highly non-uniform, and exhibit wave pat- 
terns in the x — z plane that are almost stationary over 
time. Results from 3D simulations show very similar 
features in the x — z plane. In particular, in runs R30Z1- 
3D and R10Z1-3D, there is a clear segregation of parti- 
cles with different stopping times, and their wave pat- 
terns have a phase shift relative to each other. How- 
ever, in the y — z plane, there is no coherent structure in 
the projected distribution of particles in any of our 3D 
simulations. This contrasts with the expectations from 
the KHI , where the particle layer kinks and br eaks into 
clumps (|Johansen et al.ll2006at lBarrancoll2009t ). Based 
on this observation, we infer that in our 3D simulations, 
KHI is not present in the azimuthal direction. Moreover, 
in the x — y plane, we see azimuthally elongated stripes 
of the large particles (in black). This feature, together 
with the standing wave structure in the x — z plane, is 
most likely to be due to SI. KHI resulting from the verti- 
cal shear in the gas radial velocity is another possibility, 
however, we have found that Ri x is always larger than 
Riy from our simulations, therefore the KHI is unlikely 
to play a role in the simulations presented here. 

In Figure [3] we show the Richardson number profile 
associated from u y calculated from the saturated states 
of all our simulations. The Richardson number is gen- 
erally smallest in the disk midplane, and increases with 
height. In almost all our 3D simulations (dashed curves), 
Ri y is greater than the critical value (0.1), therefore, the 
dusty midplane layer is expected to be stable against 
vertical shear, consistent with the spatial distribution of 
particles discussed above. Given the fact that Ri does 



TABLE 2 
Vertical diffusion coefficient 



Run 


Z 


A,,. 


(2D) 


O s ,*(3D) 




0.01 


2.05 x 


10~ 5 


1.51 x 10~ 5 


R41 


0.02 


1.72 x 


10" 5 


- 




0.03 


5.09 x 


10~ 6 


0.90 x 10~ 5 




0.01 


2.12 x 


10~ 5 


1.42 x 10" 5 


R21 


0.02 


1.05 x 


10" 5 


- 




0.03 


1.97 x 


10" 6 


1.57 x 10" 5 




0.01 


5.14 x 


10~ 5 


4.82 x 10" 5 


R30 


0.02 


1.14 x 


KT 4 


- 




0.03 


6.28 x 


10~ 6 


1.21 x 10" s 




0.01 


1.91 x 


10~ 4 


1.10 x 10" 4 


RIO 


0.02 


2.44 x 


10~ 4 


- 




0.03 


5.63 x 


10~ 5 


3.03 x 10~ 5 



The diffusion coefficients are measured in unit of c s H g . 

not solely determine stability, this observation does not 
entirely exclude the possibility that Ri y could be main- 
tained by KHI. However, it is important to note that 
KHI is suppressed in 2D. We see that Ri y from all our 
2D simulations (solid curves) are generally close to their 
3D counterpart. This means that the SI itself is able to 
maintain Ri above the critical value, and suggests that 
the KHI is indeed absent in all our simulations. 

The main reason that we do not observe KHI is that 
the turbulence generated from the SI is strong enough 
to prevent particles from settling sufficiently to trigger 
KHI. We note that the strength of the SI turbulence de- 
creases as the particle stopping tim e r s decreases (as ex- 
pecte d from the linear analysis of lYoudin fe Goodman! 
120051 and as confirmed by our numerical experiments). 
The turbulence in our simulations is mainly generated 
from relatively large particles with t s > 0.01 (see also 
the next subsection). We have not explored the regime 
where all particles are strongly coupled to the gas. How- 
ever, in this regime, we expect the SI to be generated on 
much smaller spatial scales with much lower amplitude, 
so that the particles settle until the KHI is triggered. In 
this regime, the dust-gas system behaves as a single fluid, 
where the dust contributes to the mass density but not 
the pressure o f the fluid. This is the approach adopted by 
iGhianel (l2008f l . iBarrancol (f2009h and lLee et all mnM O} 
to study the KHI. 

3.2. Density Profile and Vertical Transport 

Figure [4] shows the vertical density profiles for parti- 
cles of different types from all our 3D simulations, calcu- 
lated by binning the particles into vertical grid cells and 
averaging over time after saturation. Results from 2D 
simulations are generally similar. 

The vertical density profile of particles is determined 
by the balance between particle settling and turbulent 
diffusion. Unlike studies of passive particles under the in- 
fluence of homogeneous externa l turbulence (jCuzzi et al.1 
ll993HYoudin fc Lithwickll2007l ). the turbulence from our 
simulations is self-generated, and is non-homogeneous 
(strongest at the disk midplane). To study the prop- 
erties of turbulent diffusion, one approach would be to 
assume some functional form for the vertical profile of the 
diffusion coefficient D„ iZ (z), and fit the particle density 
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Fig. 3. — Profiles of the Richardson number computed from the gas azimuthal velocity from all simulations. Each panel shows the results 
from one series of runs. Red, blue and black curves label 2D simulations with Z = 0.01,0.02 and 0.03 respectively. The corresponding 
3D simulations results are mar ked with the sa me colors but using dashed lines. Horizontal dash-dotted line marks the critical Richardson 
number Ri = 0.1 adopted from [Chiang ( 2 0081 ) . 
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Fig. 4. — Vertical profiles of particle density from all our 3D simulations. Profiles of different particle types arc labeled with different 
colors, with the same color scheme as in Figure [l] and particles with smaller stopping time have more extended profiles (see Tabic [T] for 
reference). Solid lines are time averaged particle density profiles from the saturated states of our runs, while dashed curves are model 
dens ity profiles by assuming a constant turbulent diffusion coefficient fits to the density profile of particles with largest stopping time. See 
EH for details. 

profiles. However, after several experiments we found it 
difficult to fit the density profile of all particle species si- 
multaneously with any simple functional form of D gtZ (z) 
6 . In fact, the wave patterns in the x — z plane shown 
in Figure [5] suggests that the classical turbulent diffusion 
scenario may be too simple. 

Instead of fitting the vertical profile of the turbulent 
diffusion coefficient in the gas, we pose the question in 
another way: What is the effective vertical diffusion co- 
efficient at the disk midplane for the particles that are 



6 Part of the reason is that the Schmidt number Sc, defined as 
the gas diffusivity divided by the particle diffusivity, is uncertain. 
In the limit r s <C 1, one expects Sc — > 1. Even in this regime, 
we find the resulting profile Dg iZ (z) is not described by any simple 
functional form that works for all our runs. 



driving the turbulence? Since we have identified the SI 
as the source of the midplane turbulence, one expects 
particles with relatively large stopping times to drive 
the turbulence both from a theoretical point of view 
(|Youdin fc Goodmanll2005D and from non -stratified sim- 
ulations of SI dJohansen fc Youdinl 12007. Bai & Stone, 
unpublished) . To address these questions more quantita- 
tively, we find the following approach particularly useful. 
We fit the horizontally averaged vertical density pro- 
file of the largest particles r s = r max in each simulation 
using the classical picture of turbulent diffusion. Since 
these particles (as well as particles with slightly smaller 
r s ) actively drive the disk turbulence, the gas turbulent 
diffusion coefficient across this particle layer can be re- 
garded as constant. Therefore, the vertical density profile 
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of thes e particles is expected to be Gaussian, with scale 
height (|Youdin fc Lithwicld 120071 ) 



H p (t s ) 



D g ,M 



Qt* 



t s t? 



(15) 



where r e = fH c ddy is the turnover time of largest ed- 
dies. The basic assumption behind this formula is 
stochastic turbulent forcing on passive particles with 
the autocorrelation function of the turbulence Pit) = 
exp (— i/t e ddy)/27r, corresponding to a Kolmogorov spec- 
trum. We do not have much knowledge of t c( jdy for SI 
turbulence, but expect it to be comparable to the or- 
bital time (the only time scale of the problem) , and take 
r e = 1. The exact value of r e does not matter much, 
since it only gives an order unity correction to H p . 

By fitting the vertical density profile of the largest par- 
ticles with a Gaussian we obtain D giZ (0) for all simula- 
tions, and the results are summarized in Table [2] For 3D 
runs, the results are also plotted in Figure 0] as dashed 
lines. We see that the vertical profiles of the largest 
particles are well fitted with a Gaussian. In addition, 
we predict the vertical density profile for other particle 
species, using equation (|15j) and assuming a diffusion co- 
efficient which is constant with height. Obviously, this 
will overpredict the scale heights for small particles, since 
they respond to the turbulence passively. However, for 
particles that actively participate in the instability, we 
expect their density profile to be comparable to the pre- 
dicted profile, since they are driving turbulence to main- 
tain D gz close to Dg^ z (0) across their scale heights. In 
this way, we are able to identify the particle species that 
are responsible for the disk turbulence (hereafter termed 
as "active" particles). 

From the R41 runs, we see that active particles range 
from t s = 0.1 (for R41Z1) to r s > 0.01 (R41Z3). Active 
particles for R21 runs have r s > 0.03. For R30 runs, 
particles with r s > 0.03 are active, while for R10 runs, 
all particles are active. We see that although there is 
a diversity in the size range of active particles, which 
depends on both solid abundance and particle size dis- 
tribution, the minimum size of active particles for most 
of our runs is about t s = 0.01 — 0.03. For run R41Z1, 
although we have identified somewhat larger r s values 
for active particles, particles with t s = 0.01 — 0.1 must 
actively participate in the instability because the abun- 
dance of t s =0.1 particles alone is too small to trigger 
SI. 

Next we study the midplane diffusion coefficient from 
our simulations. We emphasize that the strength of the 
turbulence (hence D 9yZ ) is self-regulated: the settling of 
particles continues until the turbulence they generate is 
sufficient to stop the settling. To better interpret our 
results, we construct a toy model describing the self- 
regulated turbulence. In this model, we assume all par- 
ticles are active, and that the particles are single-sized, 
with fixed stopping time r s . Since all particles are active, 
their vertical density profile can be approximated by a 
Gaussian, so that the particle to gas mass ratio at the 
disk midplane is given by e = ZH g /H p . The midplane 
diffusion coefficient D„ :Z depends on both r s and e. For 



simplicity, we parameterize the dependence as 

D g , z = e a f(r s )H 2 g n , (16) 

where D is normalized to H g tt, /(r s ) is a coefficient that 
incorporates the dependence of D on r s , and a is a power 
law index that reflects the sensitivity of the dependence 
of D on e. We note that D gz — > at both e — > and 
e — > oo, therefore, we expect a > when e is small and 
a < for large e. Using equation (fi"5j) and neglecting the 
second square root (which is order unity) on the right 
hand side, we obtain 



J Z a ) H, 



H„ 



D g , z =(fZ Q )^T. 
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?$= ■ H 2 g n 
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In the above equations, the dependence of particle scale 
height and diffusion coefficient on Z is reflected in the 
index a. When a is positive, increasing Z leads to larger 
H p and larger D 9iZ . When a is negative, the situation re- 
verses. Below, we apply this simple model to our results. 
Since our simulations contain multiple particle species, 
we may take e to represent the contribution from all par- 
ticle species participating in the SI (i.e. with t s > 0.01). 

Our R30 and R10 runs show similar behavior between 
2D and 3D simulations with respect to vertical diffu- 
sion properties. Increasing Z from 0.01 to 0.02 produces 
stronger turbulence, while further increasing Z to 0.03 
dramatically reduces H p . This corresponds to the tran- 
sition from a>0toa<0ata threshold e (hence thresh- 
old Z = Z t h). Beyond Z t h, H p sensitively depends on Z 
because the corresponding power law index a/(2 + a) 
quickly drops to large negative values once a turns nega- 
tive. Consequently, a small increase in Z results in strong 
particle settling and greatly enhances midplane particle 
density. This result has important implications for par- 
ticle clumping discussed in the next section. 

In our 2D R41 and R21 runs, we see that D gz mono- 
tonically decreases with Z, suggesting a < for Z > 
0.01. Based on this result, we infer that the strength of 
the SI for a particle size range r s = 0.01 — 0.1 is a de- 
creasing function of e for e > 0.5. The 3D simulations 
give somewhat different results. For both 3D R41 and 
R21 runs, D gz slightly increases with Z at least in the 
range Z < 0.03, indicating a > 0. It is very likely that 
the threshold abundance Z t h is above 0.03, which is sub- 
stantially larger than their 2D counterparts. We note 
that the behavior of the SI turbulence for t s < 0.1 par- 
ticles in 3D is different from that in 2 D in non-stratified 
simulations (jJohansen fc Youdinll2007n . Our results indi- 
cate that the difference remains when vertical gravity is 
included, and 3D simulations are needed to better catch 
the dynamics of small particles. 

In our toy model, all of our ignorance on the depen- 
dence of -Dg, z on t s is encapsulated in the unknown func- 
tion f(r s ). From Table H we sec that the R30 and R10 
runs generally have larger D gz than R41 and R21 runs. 
This result implies that turbulence generated from larger 
particles r s ~ 1 is stronger than that from smaller par- 
ticles, i.e., /(t s ) is an increasing function of t s in this 
range, consistent with results fro m non-stratified simula- 
tions (jJohansen fc Youdinl [20071 ). 
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In sum, we have identified that particles actively par- 
ticipating in SI generally have stopping time r s > 0.01. 
The strength of the turbulence largely depends on the 
density of these active particles at disk midplanc. We 
find that the particle scale height (thus the turbulent dif- 
fusion coefficient) strongly depends on solid abundance. 
Such strong dependence is caused by a sharp drop in the 
strength of the turbulence with increasing particle to gas 
mass ratio e when e is larger than a certain threshold 
value. 

4. PARTICLE CONCENTRATION 

4.1. Formation of Particle Clumps 

Probably the most interesting property of the SI is the 
concentration of particles. The degree of particle con- 
centration strongly depends on the mass distribution of 
solids in PPDs. In our simulations, we normalize particle 
density to the background gas density at the disk mid- 
plane pg t b{r,z = 0). A useful scale to measure particle 
concentration is the Roche density, above which the par- 
ticle clump can be consid ered as gravitationally bound 
(jBinnev fc Tremaind[2008l) 



Piochey) 



3M* 



1.34x10 



3 (/ a/ /t) 1/2 2» 



./•, 



r AU 2 PgA r ,°)- 

(18) 

The normalized Roche density (relative to the back- 
ground gas density at midplane) scales as the square root 
of stellar mass and disk temperature, and is inversely 
proportional to disk mass, meaning that the Roche den- 
sity is easier to reach for massive disks (with large f g ). 
In the MMSN model, the Roche density is of the order 
proche = 10 3 £o,b: an d only weakly depends on r as r -1 ' 4 . 

In Figure [5] we show the time evolution of maximum 
particle density p Pjmax from all our simulations. We 
first look at results from 2D simulations. For all the 
four run series, /Op ima x increases with solid abundance 
Z. However, the dependence of p Pimax on Z is highly 
non-linear. For run series R21, R30 and RIO, there is 
no significant clumping of particles for Z = 0.01 and 
0.02. However, significant clumping occurs at Z = 0.03, 
with maximum particle density reaching 10 3 times the 
background gas density, comparable to the Roche den- 
sity (IT51). This trend is consistent with the results by 
iJohansen et al.l (120091) (see also the supplemental infor- 
mation in IJohansen et al.l 120071 ) , who considered parti- 
cles with stopping time in the range of r s = 0.1 — 0.4. 
As emphasized in the previous section, there is a sharp 
enhancement of averaged midplanc particle density with 
increasing Z once Z exceeds some threshold value. This 
density enhancement further favors strong concentration 
of particles by SI, which explains the trend we have ob- 
served in Figure [5J 

The particle clumping also depends on the particle size 
distribution. In the R41 run series, where the majority 
of the particle mass resides in strongly coupled particles 
t s < 10~ 2 , we see that there is no significant clumping of 
particles up to Z = 0.03. As noted in the previous sec- 
tion, particles that effectively participate in SI are those 
with relatively large stopping times t s > 10~ 2 . These 
particles are also the ones that actively participate in 
the clumping (see the next subsection). For R41 runs, 



the abundance of these "active" particles is much smaller 
than our R21, R30 and R10 runs, which makes the criti- 
cal (total) abundance for strong particle clumping larger. 
In fact, we do observe strong clumping when we increase 
the total abundance to Z = 0.05. Based on the dis- 
cussion above, we conclude that in order for the SI to 
efficiently concentrate particles, the mass of the solids 
with stopping time t s > 10~ 2 should exceed a critical 
value Z cv - lt . The results from 2D simulations suggest that 
S T >io- 2 ^ k ~ -^crit ~ 0.02 is necessary for significant 
particle clumping 7 . 

The 3D simulations show similar trends as in 2D, 
but the condition for strong particle clumping is more 
stringent. Among the eight 3D runs, strong clumping 
occurs only in run R10Z3-3D. The maximum density 
for all other runs remain small in the saturated state 
(pp,max ^ 50/Og.f,). In particular, the 3D R21Z3 and 
R30Z3 runs do not show clumping as in their 2D counter- 
parts, and both of them have larger D gz . Since KHI is 
unlikely to be present in these simulations, the different 
results between our 2D and 3D simulations should be at- 
tributed to the different behavior of the SI in 2D and 3D. 
It appears that the formation of dense particle clumps fa- 
vors the mass distribution of particles to be dominated 
by larger particles than in 2D, or larger values of Z CT n is 
needed. 

Interestingly, in run R30Z3-3D, a very dense clump 
(actually a nearly axisymmetric stripe) forms at about 
t = 150f2 -1 . The composition of this (transient) clump 
is similar to its counterpart R30Z3-2D (see next subsec- 
tion). It lasts for about 10 orbital times and then is grad- 
ually dissolved. Both the Richardson number profile and 
particle distribution disfavor the presence of KHI during 
the process. Nor is there any significant vorticity gener- 
ation in the vicinity of the clump which might indicate 
KHI. By comparing with Figurc[T] we see that the period 
during which the clump is dissolved is accompanied by 
an increase of the height of relatively small particles with 
t s < 0.1. It is likely that the formation of the transient 
clump is due to our unrealistic initial condition 8 . 

The results we have obtained show a clear dichotomy 
on the particle concentration properties. Specifically, 
the maximum density is either very small with p v>max < 
50p g .b, or very large with p P!max > 1000p Sjb . Self-gravity 
becomes important when the particle density approaches 
the Roche density (fT5)l . This means that for our simu- 
lations that do not show signature of strong clumping, 
adding self-gravity will not change the picture qualita- 
tively 9 . For simulations with strong clumping, the max- 
imum particle density is already comparable with the 

7 The value of the critica l metallicity also dep ends on the pres- 
sure gradient parameter 14 (Bai & Stone 2010b). A smaller value 
of El leads to smaller Z cr it ■ 

8 As small particles diffuse towards larger heights, the gas az- 
imuthal velocity at disk midplane is reduced, thus larger particles 
feel a stronger headwind, enhancing the turbulence strength of the 
SI, which destroys the clumps. 

9 Recent N-body simulations by Michikoshi et al. (2010) show 
that gravitational collapse may occur before Roche density is 
reached due to the drag force. This is unlikely to affect our con- 
clusion because in the non-clumping case p p , m ax is usually more 
than one order of magnitude smaller than the Roche density, and 
densest regions are only transient. 
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Fig. 5. — The evolution of niaxinium particle density for all our simulations. Each panel shows the results from one run series. In the 
upper panels, red, blue and black curves label 2D simulations with Z = 0.01, 0.02 and 0.03 respectively. The 3D results are shown in the 
lower panels with Z = 0.01 and 0.03 marked with red and black. The maximum density is normalized with respect to the background gas 
density at the disk midplanc. 

Roche density, and in this case we expect the forma- 
tion of a few p l anete simals from the simulations as in 
Uohansenet all (|2009f ). 

Particle concentration properties are known to depend 
on numerical resolution. To assess the validity of our re- 
sults, we have also performed the same set of simulations 
with half our standard resolution. We find the same di- 
chotomy between strong clumping and no clumping. The 
only exception is the R30Z3-3D run: it shows strong par- 
ticle clumping in the low-resolution run which does NOT 
dissolve as in our standard resolution run. The reason 
is that the turbulence generated from the lower resolu- 
tion run is weaker, thus particles settle more which favors 
clumping. This test justifies the necessity of conducting 
high resolution simulations. In the mean time, it sug- 
gests that the critical abundance for particle clumping 
in this run may be only slightly larger than 0.03. There- 
fore, the particle clumping properties from 2D and 3D 
simulations is not dramatically different when r max = 1. 

4.2. Properties of Dense Clumps 

In this subsection we discuss more details of the three 
simulations that exhibit strong particle clumping. First, 
we examine the composition of these dense clumps by 
plotting the cumulative probability distribution function 
(CPDF) of particle densities for different particle species 



P(p p > p). The CPDF measures the probability of a par- 
ticle residing in a region with total particle density larger 
than p. In Figure[6j we plot the CPDFs of the three runs: 
R21Z3, R30Z3 and R10Z3. At relatively high densities 
with p p > 10 2 p g , we see that in all three cases, the dense 
regions are composed of particles with the largest stop- 
ping times. In run R21Z3-2D, the mass fraction of differ- 
ent particle species in the dense clumps is increasing with 
particle stopping time t s , and is completely dominated 
by the largest particles t s = 0.1. In the case of R30Z3-2D 
and R10Z3-2D, where the largest particles have r s = 1.0, 
the composition of the clumps are dominated by the two 
largest particle species. Contribution from other particle 
species to the clumps is almost negligible by mass. 

For R21Z3 and R30Z3 runs, 3D simulations do not 
show particle clumping, therefore, the resulting CPDFs 
differ substantially from those in 2D runs. Nevertheless, 
these CPDFs provide typical examples for simulations 
without clumping. The shapes of the CPDFs from differ- 
ent particles are very similar, and curves for larger parti- 
cles are located to the right of those for smaller particles, 
consistent with the vertical stratification of particles. For 
run R10Z3-3D, the particle clumping is stronger than the 
2D case, and the densest clumps are almost equally made 
of particles with r s = 1 and t s = 10 -1 / 2 . 

Next, we consider the motion of the dense clumps. In 
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Fig. 6. — The cumulative probability distribution function (CPDF) of particle densities in the saturated state from the three simulations 
that exhibit strong particle clumping. Different colors mark the CPDF of different particle species, as labeled in the legend. Solid curves 
come from 2D simulations, while results from 3D simulations are denoted in dashed curves. Particle densities are normalized to background 
gas density at disk midplane. 




Z=0.01 



A t=60tl' 




A t=9oa 



-1 





-20 -15 -10 -5 
xftr) 



-40 -30 -20 -10 -60 -40 -20 

x ft r) x (r, r) 



10" 



c 

CD 
T3 



>.10 



.Q 
CO 
.Q 
O 



a. _ 2 
10 



A t=3on _1 


1 


Z=0.03 

A t=60ri" 1 


I 




_T=1(T 1 

x=10-°' 5 
— — s 

T=1 


At=90£2 1 




- 




\ 












/ 


^ 


V( 




















/ 


If 


I 






V 


rl , J f 


J 






I. 


I 


V 


l] 


J 


r 



-15 



-10 -5 

x(tir) 



-30 -20 -10 -50 -40 -30 -20 -10 

x ft r) x (ti r) 



Fig. 7. — The probability distribution of the radial distance traveled by various types of particles after a given time interval AT from 
our runs R10Z1-3D (top) and R10Z3-3D (bottom). We choose At = 30, 60 and 90 as shown in left, middle and right panels respectively 
(unit is f2 _1 ). Particles of different types are labeled by different colors (see the legend). 



FigureEl we mark the location of the densest point with a 
red dot in runs with strong particle clumping. By moni- 
toring the location of the densest point with time, we find 
that it wanders slowly. Another useful way of studying 
the dynamics of the clumps is by tracking the radial tra- 
jectories Xi(t) of a sample of particles. We relocate the 
particle positions when they cross the radial boundaries 
of our simulation box so that their trajectories are con- 
tinuous. By tracing a large number of particles in the 



saturated state of our runs, we obtain the distribution of 
x(t + At) — x(t) for each particle species at time interval 
At. In Figure [7] we show the probability distribution of 
x(t+At) — x(t) for a number particle species from our run 
R10-3D. When Z = 0.01, no particle clumping occurs. 
The distribution of x(t + At) — x(t) is close to a Gaus- 
sian (or a parabola in logarithmic scale) and the width 
increases with At, consistent with undergoing a random 
walk. Meanwhile, the center of the distribution drifts in- 
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ward with time (see ^Sl for more discussion). However, 
when particle clumps are present, as in the Z = 0.03 
case, the shape of the distribution deviates substantially 
from a Gaussian, especially for particles that make up the 
clumps (the largest particles, shown in the blue and green 
curves). For these clump-making particles, the width of 
particle distribution still increases with At, as expected 
from turbulent diffusion, but a substantial fraction of 
theses particles stay nearly stationary without drifting 
(near x = 0), making the resulting distribution more 
and more elongated with time. The leftmost location of 
the particle distribution moves inward with time, and is 
set by the radial drift velocity. More interestingly, we 
see almost evenly separated multiple peaks in the distri- 
bution function. In fact, the separation between these 
peaks equals the radial size of our simulation box. The 
physical picture becomes clear that the clumps stop some 
of the particles from drifting radially, and particles are 
kept in the clump for a few orbits or more before leav- 
ing for the next clump. Similar behavior is observed for 
other runs with particle clumping. 

5. RADIAL TRANSPORT OF SOLID PARTICLES 

As expected from particle-gas equilibrium, particles ex- 
perience head wind from the gas and drift radially in- 
ward. Particles with different stopping times drift at dif- 
ferent velocities. At the same time, the instabilities gen- 
erated at the disk midplanc diffuse the particles. These 
two processes transport particles radially in PPDs, and 
is the subject of this section. In particular, we show that 
it is important to study the radial transport of particles 
by considering particles of all sizes simultaneously, rather 
than individually. 

5.1. Radial Drift Velocity 

We calculate the averaged radial drift velocities for 
each particle species from all our runs, and the results 
are shown in Figure |S] The measured mean drift ve- 
locities are shown in squares (2D) and circles (3D). We 
have also plotted the la limits for particle drift veloc- 
ity based on the rms fluctuations, which are indicated in 
blue and red vertical bars. In the figure, the velocities 
are normalized to r/vx. Clearly, the radial drift velocity 
monotonically decreases with particle stopping time, and 
the drift is fastest for marginally coupled particles. 

The classical result on the radial drift of particle s is 
the NSH equilibrium solution (iNakagawa et al.lll986f ). It 
describes the equilibrium state between solids and gas in 
unstratified (neglecting vertical gravity) Keplerian disks, 
where gas is partially supported by radial pressure gra- 
dient. In the NSH equilibrium, the drift speed is given 

by 

2r s 



(1 + ef + t\ 



;V V K 



(19) 



We emphasize that the conventional NSH solution is ob- 
tained by considering a single species of solids. Equation 
HU) does not simply generalize to the case with multiple- 
species of particles by replacing e to e^ for each particle 
species k. In Appendix [A] we provide the generalized 
formula for multi-species NSH equilibrium, and the so- 
lution involves evaluation of an inverse matrix of order 



2iVt ypc . It reflects the fact that although different parti- 
cle species do not interact directly with each other, they 
are indirectly coupled via their interactions with gas. 

In Figure [3 the bold solid lines show the expected ra- 
dial drift velocities from single-species NSH equilibrium. 
We see that there are large deviations from the measured 
mean drift velocities, with two notable features. First, 
for relatively large particles, the drift velocities are re- 
duced from single-species NSH values. The reduction 
is strongest for runs with the largest Z. Second, the 
smallest particles drift outward, rather than inward as 
expected from the single-species NSH solution. 

To calculate the expected radial drift velocity from a 
multi-species equilibrium, we first use the particle den- 
sity profiles extracted from £ )3.2I and calculate the drift 
velocity in each vertical bin. The drift velocity is then 
weighted by particle density in each bin to yield the 
mean drift velocity. The results are plotted in dashed 
and dash-dotted lines (for 2D and 3D runs respectively) 
in Figure \8\ We see that these curves provide an excel- 
lent fit to the measured mean radial drift velocities in all 
simulations. In fact, the two features mentioned above 
are natural consequences of the multi-species solution. 
Due to the sub-Kcplcrian motion of the gas, particle 
drag increases gas angular momentum, leading to out- 
ward drift of gas. In the presence of both weakly coupled 
and strongly coupled particles, the strongly coupled par- 
ticles are tied to the gas and therefore drift outward with 
the gas. Marginally coupled particles still drift inward, 
but due to the influence of the smaller particles, these 
particles feel a weaker headwind (i.e., the gas azimuthal 
velocity is closer to the Keplerian value), resulting in a 
smaller drift velocity compared with the single-species 
solution. With increasing Z, thus higher midplane parti- 
cle density, the gas becomes more entrained by the solids, 
leading to stronger reduction of the drift velocity for large 
particles. 

The residuals from the multi-species NSH solution fit 
to the measured mean drift velocities are largest for par- 
ticles with largest r s , likely due to their participation 
in SI, and/or clumping . In th e non-stratified simulation 
of iJohansen fc Youdinl (|2007t ) , it was shown that in the 
saturated state of SI, the radial drift velocity is either in- 
creased or decreased depending on run parameters. In 
our simulations, these effects are secondary compared 
with the multi-species effect. The measured drift ve- 
locities from 2D (squares) and 3D (circles) simulations 
generally agree with each other. The (small) differences 
can be attributed to the differences in the particle verti- 
cal density profiles. 

So far we have focused on the mean radial drift ve- 
locities. In the saturated state of our simulations, the 
particle radial drift velocities follow a distribution, due 
to the SI. We see in Figure [5] that in most of the runs, 
the fluctuation level is about (0.05 — 0.15)t]Vk ■ This fact 
is closely related to the radial diffusion of particles dis- 
cussed in the next subsection. Based this observation, 
we can estimate the particle radial diffusion coefficient 
to be D x ~ {0.1r]v K ) 2 /n ~ 2.5 x lQ- 5 c s H g . 



5.2. Radial Diffusion 
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Fig. 8. — Radial drift velocities for different particle species from all our simulations. Blue squares: averaged particle radial drift velocity 
from 2D simulations, and vertical blue bars indicate its 1<t limits. Red circles: averaged particle radial drift velocity from 3D simulations, 
and vertical red curves show its Iff limits. We have slightly shifted the symbols for 2D and 3D runs in order to show the la bars more 
clearly. Bold black solid line: radial drift velocity expected from single-species NSH equilibrium. Black dashed (dash-dotted) line: radial 
drift velocity expected from multi-species NSH equilibrium (see Appendix [XJ for 2D (3D) simulations. All velocities are normalized to 
rjv K . 

we discuss these results further. 

First, the above procedure for measuring the diffusion 
coefficient does not apply to runs that show strong par- 
ticle clumping. As we see in Figure [3 the distribution of 
x(t + At) — x(t) deviates strongly from a Gaussian due 
to the influence of the clumps. The measured width of 
the distribution is about half the distance traveled by the 
fastest drifting particles (those that are not confined in 
the clumps), and we observe that er 2 scales as At 2 rather 
than At from our measurement. Therefore, the measured 
D x from R21Z3-2D, R30Z3-2D and R10Z3 (both 2D and 
3D) runs for those clump making particles (or the largest 
two particle species in the run) is not valid. In Figure HI 
we see the measured D x for these particles have anoma- 
lously large values. Such particles can reside in the disk 
for much longer than if there were no clumping. 

Next, we discuss diffusion of non-clumping particles. 
In each simulation the measured D x generally approaches 
an asymptotic value for particles with r s < 10 -2 , but is 
different between different particle species for particles 



The radial diffusion of particles is generally character- 
ized by the radial diffusion coefficient D x . From our sim- 
ulations, we can measure D x for different particle species 
based on the random walk model of particle diffusion. 
We calculate the distribution of shift in the particle ra- 
dial position at various time intervals At as in Figure 
[7J and measure the width (rms) of the distribution a as 
a function of At. The spreading due to a random walk 
results in an Gaussian distribution, and a is related to 
the diffusion coefficient by 



x ~ 2 dt ' 



(20) 



For each particle species, we measure a 2 for different 
At, and fit the slope in of the a 2 — At curve by linear 
regression. The results are summarized in Figure [9] The 
range of the radial diffusion coefficient is consistent to 
within an order of magnitude of the estimate in the last 
subsection based on the spread of radial drift velocities. 
It is also comparable with the vertical diffusion coefficient 
at disk midplane estimated in £13.21 (sec Table [2]). Below 
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Fig. 9. — Radial diffusion coefficient for different particle species from all our simulations. Results from 2D and 3D simulations arc shown 
in solid and dashed lines respectively. Red, blue and black curves represent different metallicities with Z = 0.01, 0.02 and 0.03 respectively. 
Diffusion coefficients are normalized by c s H g . 



with r s > 10~ 2 . This can be due to multiple reasons. 
First, similar to the vertical diffusion of particles, the 
radial diffusion coefficient also depends on the vertical 
position in the disk, and the radial diffusion in the disk 
midplane is expected to be the strongest. Our measured 
D x can be considered as a vertically averaged quantity. 
Therefore, D x is expected to be larger for particles with 
larger t s , since they stay closer to the midplane. This 
trend is observed in runs R41 and R21. Second, different 
particles react differently to the turbulence. In the case 
of Kolmogorov turbulence, the particle d iffusivity scales 
as (1 + t s 2 )- x (jYoudin fc Lithwickl 120071) . This may be 
responsible for the decrease of D x towards r s = 1 in 
R30 and R10 runs with Z = 0.01 and 0.02. Thirdly, 
different particles participate in the SI in different ways 
(i.e., actively or passively). The SI may strongly affect 
the transport properties of the active particles, with the 
extreme example being the clump-making particles dis- 
cussed above. Despite the different values of D x for dif- 
ferent particle species, one may take the asymptotic value 
of D x as measured from the smallest particles as charac- 
teristic of the radial diffusion coefficient in the gas. These 
asymptotic values correlates with the vertical diffusion 
coefficient well (see Table [2]). 



To address the effectiveness of radial diffusion com- 
pared with radial drift, we denote the mean radial drift 
velocity to be v r = kt\vk , and the diffusion coefficient to 
be D x = (f3r]v k) 2 / Q ■ After time t, the ratio 



c,.f 



'-^-f^- 



(21) 



reflects the relative importance between radial drift and 
turbulent diffusion, where a{t) = ^/2D x t. Diffusion is 
important when £ < 1 . From Figure |8l we see that for 
the largest particles, k > 0.1. From Figure O we have 
/? < 1. Therefore, the effect of radial diffusion of par- 
ticles becomes negligible compared with radial drift be- 
yond 100 orbital periods. Again, this discussion does not 
apply to the situation when particle clumping is present, 
where large particles can be retained in the clumps and 
some of them may survive the radial drift. 

6. COLLISION VELOCITIES 

The initial stage for planetesimal formation is the 
growth of solid bodies by mutual collisions. The size 
distribution of particles in the PPDs therefore depends 
on the outcome of two-body collisions, which further de- 
pends on the properties of the colliding particles (e.g., 
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size and porosity) and collision velocity. Laboratory 
experiments show that at low collision velocities (< 
lm-s -1 ), collisions generally lead to sticking or bounc- 
ing. Larger collision veloc ities tend to result in fragmen- 
tation (see the review bv lBlum fe Wurmll2008l ). Never- 
theless, sticking can also occur with collision velocities 
up to 10 — 20m-s~ 1 in some regimes (see Figure 11 of 
iGiittler et al.l|2009| ). The particle size distributions used 
in this paper can be considered as a first approximation 
to the outcome of grain growth in PPDs. In turn, we 
can measure the two-body collision velocity produced by 
the SI from our simulations and investigate whether our 
selected particle size distribution is consistent with the 
outcome of collisional coagulation. 

We measure the relative speeds of all particle pairs 
within a distance Ar in the saturated state of our simu- 
lation snapshots. These velocities form a representative 
sample of particle relative velocity distribution (RVD) in 
the vicinity of a tracer particle. We assume that par- 
ticles that collide with this tracer particle would have 
the same RVD. The measured RVD depends somewhat 
on the choice of Ar. In practice, we choose Ar to be 
a quarter of a cell size, in order to reduce the (misrep- 
resented) measured collision velocity between strongly 
coupled particles (see Figure [10] and the discussion that 
follows), while maintaining good statistics. To obtain the 
distribution of collision velocities with a tracer particle, 
the RVD must be weighted by the relative velocity, since 
the collision frequency is enhanced at larger relative ve- 
locities. The corresponding CPDFs (similar to £|4.2[) arc 
shown and discussed in Appendix [Bj In this context, 
it measures the probability of a particle that undergoes 
collision with relative velocity greater than a given value. 
Particle velocities are normalized to the gas sound speed 
c s in our simulations. In all the results presented in this 
section, we adopt c s = 0.99km s" 1 , corresponding to the 
MMSN model at 1AU. 

In order to visualize the particle collision velocities in 
a compact way, we characterize the CPDFs by the me- 
dian collision velocity (at P = 0.5) and its lc limits (at 
P = 0.68 and P = 0.32). In Figure QUI we show the me- 
dian collision velocities and ltr limits for various pairs of 
particle species from all our 3D simulations. Results from 
2D simulations are generally similar, and are not plotted. 
To interpret these results, we consider two sources of the 
collision velocities: radial drift and turbulence. 

To calculate the contribution from radial drift, we eval- 
uate the multi-species NSH equilibrium in each vertical 
cell bin j (j — 1 , . . . , N z ) , from which we obtain the rel- 
ative radial drift velocity (Au r )j. k between each pair 
of particle types k\ , hi in that bin. The relative veloc- 
ity is further weighted by collision frequency in that bin, 
proportional to (Av r ) k k e k e k . Integrating over all the 
vertical bins, we obtain the expected collision velocity 
from radial drift, which is shown as solid curves in Fig- 
urellOl We see that with the exception of run R10Z3-3D, 
these curves fit the median collision velocities very well, 
meaning that relative radial drift is the dominant source 
of collision velocities. 

R10Z3-3D is the only 3D run that shows strong parti- 
cle clumping, and the measured median collision velocity 



is strongly reduced from our predictions. This is clearly 
seen in the CPDF plot (see Figure [12] in Appendix [E| . 
However, in these simulations, the median collision veloc- 
ity no longer characterizes the overall collision velocities 
because the shapes of the CPDFs are strongly deformed 
due to the clumping. In fact, there is still a high- velocity 
tail in the CPDF of collision velocity, which reaches val- 
ues as high as 30m s _1 . This tail is most likely caused 
by collisions outside the clump, as indicated in Figure 
1131 and our predicted collision velocities should apply in 
these low density regions. 

The relative radial drift velocity can not account 
for the collision velocity between particles with the 
same stopping time (therefore all solid curves reach a 
zero point in Figure 1 10[) . To remedy this limitation, 
we further consider the contribution from turbulence. 
So far turbulence induced particle collision velocities 
has been studied theoretically only in the framework 
of passive particles in uniform Kolmogor ov turbulence 
()Voelk et al.lll980HMarkiewicz et aT1ll991[ ). and in MRI 
turbulence (]Carballido et al.l 120081 ). We consider the 
clo sed form expres s ion o f turbulent collision velocities 
by lOrmel fc Cuzzil (|2007f ). which is based on the Kol- 
mogorov spectrum. Although these assumptions do not 
quite apply in our simulations, we adopt this approach 
as an approximate treatment of turbulence induced col- 
lision velocities. We use their equation (16), and more 
specifically, we fix the turn over time for the smallest 
eddy to be t n = 0, and take y* = t*/t s top = 1-6 as an 
approximation (where t* is the turn over time of the criti- 
cal eddy with which the particle in question is marginally 
coupled). The turn over time of the largest eddy t^, is 
considered as a fitting parameter 10 . Because the strength 
of the turbulence is vertically stratified, we take the av- 
eraged radial diffusion coefficient D x from the smallest 
particles in each of our simulation run. The averaged gas 
velocity V g is then related to D x by D x ~ VgtL- 

In Figure [TO] we also show the contribution from tur- 
bulence induced relative velocities as dashed curves. In 
order to fit the collision velocity for pairs of large parti- 
cles t s > 0.1, we find flt.L — 2 — 3 for R41 and R21 runs, 
and QtL ~ 4 for R30 and R10 runs. With this contri- 
bution, the collision velocity between the same types of 
particles can be fit very well, and it also improves the 
fit to collision velocities between particles with different 
types. 

In our R41 and R31 runs, the predicted collision ve- 
locities almost reach zero for collisions between parti- 
cles with t s < 10~ 3 , since contributions from both radial 
drift and turbulence rapidly decrease with stopping time. 
The measured collision velocities are always larger than 
the predicted values, as seen in the leftmost four panels 
of Figure [TO] and decrease towards a small asymptotic 
value at smallest r s . We have experimented with choos- 
ing different Ar in our calculations and found that the 
asymptotic value roughly scales linearly with Ar when 
Ar is less than grid size, because the gas velocity is not 

10 In principle, t^ is the same as i e< jdy defined in equation 1151 . 
where the latter is set to f2 for simplicity. Given the large uncer- 
tainties in this rough treatment of the turbulence induced collision 
velocity calculation, we allow t^ to vary. 
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Fig. 10. — The median collision velocity for two-body particle collisions from all our 3D simulations. In each subplot, upper panels 
show results from metallicity Z = 0.01 and bottom panels for Z = 0.03. Plotted in each panel are median collision velocities between one 
particle with fixed r s and a second particle as a function of r s of the second particle. Squares: measured median collision velocities from 
the simulations. Vertical bars: the la range of the collision velocities. Solid lines: expected collision velocity calculated from the radial 
drift velocities using multi-species NSH equilibrium (see Figure [8). Dashed lines: expected collision velocity from both radial drift and 
turbulence based on results from Figures \E\ and [9] (sec text for details). The collision velocities scales linearly with the adopted sound speed, 
and we take c s = 0.99km s — 1 appropriate for MMSN model at 1AU. 



resolved at scales less than a grid cell. 

From Figure I1Q[ the median collision velocity is typ- 
ically a fraction of tjvk (~ 50m s _1 with our chosen 
scaling) . Since the collision velocity is dominated by the 
radial drift, and the radial drift is largest for marginally 
coupled particles with r s ~ 1, we see that the collision ve- 
locity is relatively small in the R41 and R21 runs (where 
T max = 0.1), typically smaller than Q.lrjVK. The col- 
lision velocities from the R30 and R10 runs are much 
higher. Moreover, by comparing runs with the same par- 
ticle size distribution but different solid abundance, we 
see that the collision velocity is reduced at larger Z. This 
is again due to the reduction of radial drift velocity at 
larger Z (see Figure [8]) . The typical value of the colli- 
sion velocity in our Z = 0.03 runs are within 3m s _1 for 
R41 and R21 runs, and within 12m s" 1 for R 3 and R10 
runs. Looking at Figure 11 of iGiittler et al.l (|2009f ). al- 
though collisions with relative velocity above lm s~ x arc 
destructive in a number of situations, in other cases (e.g., 
when a porus particle hits a compact particle), particle 
growth is still possible by mass transfer with collision 
velocities less than 10 — 20m s" 1 . Detailed modeling of 
particle size evolution is beyond the scope of this paper. 
Based on the results shown in Figure [Til it is possible 
for particle growth in all our R41 and R21 runs, as well 
as R30 and R10 runs with Z = 0.03, meaning that the 
adopted particle size distribution in these runs may be re- 
alizable. On the other hand, our R30 and R10 runs with 
Z = 0.01 and Z = 0.02 appear unlikely to be realized 



in nature, due to the destructive collisions at velocities 
beyond 30m s -1 . Combined with the results in 21 we 
conclude that larger solid abundance favors grain growth 
in PPDs, which further promotes particle clumping. 

7. DISCUSSION 
7.1. Summary of Main Results 

The main purpose of this paper is to study the dy- 
namics of solids and gas in the midplane of PPDs using 
hybrid simulations. The solids and gas are coupled aero- 
dynamically, characterized by the dimensionless stopping 
time t s = flt s top- We consider a wide size distribution 
of solids as an approximation to the outcome of grain 
growth in PPDs, ranging from sub-millimeter to meter 
size. The key ingredient of our simulations is the inclu- 
sion of feedback from particles to gas. Feedback is im- 
portant when the local particle to gas mass ratio exceeds 
order unity. Moreover, it is essential for the generation 
of SI and KHI. In our simulations, we assume no external 
source of turbulence, as an approximation for the dead 
zone of PPDs. Turbulence in the disk midplane is gen- 
erated self-consistently from the SI (driven by the radial 
pressure gradient in the gas) and/or KHI (driven by ver- 
tical shear). Our simulations are local, since very high 
numerical resolution is essential to resolve the SI and 
KHI. Self-gravity is ignored, as we focus on the particle- 
gas dynamics before the formation of planctcsimals. 

Our simulations are characterized by three sets of di- 
mensionless parameters, namely the particle size distri- 



20 



Bai & Stone 



bution Tfe, solid abundance Z, and a parameter II char- 
acterizing the radial pressure gradient. In this paper, 
we fix II = 0.05, as appropriate for a wide range of 
disk model parameters (see £|2.2[) . The dependence of 
the particle clu mping properties on II is presented in a 
separate paper (jBai fc Stondl2010b| ). We consider a flat 
mass distribution in logarithmic bins in r s , and vary Z 
from 0.01 to 0.03 (see Tablc[p. We conduct both 2D and 
3D simulations, where 2D simulations are performed in 
the radial- vertical plane in order for the SI to be actively 
generated. We run the simulations for 40 — 200 orbits 
and study the properties of the particles and gas in the 
saturated state. The main results are summarized below. 

1. SI plays the dominant role in the dynamics of PPD 
midplane when the largest solids have stopping 
times t s > 10 -2 . Particles with r s > 10~ 2 ac- 
tively participate in SI, while smaller particles be- 
have passively. KHI is not observed in all our sim- 
ulations, which suggests that it may be important 
only when all particles have t s < 10~ 2 . 
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2. The strength of the turbulence generated by the 
SI and the scale height of the particle layer are 
self-regulated. There exists some threshold solid 
abundance, above which increasing Z will result 
in weaker turbulence, which promotes particle set- 
tling, leading to rapid drop of the thickness of the 
particle layer and strong particle clumping. 

3. SI can concentrate particles into dense clumps with 
solid density exceeding the Roche density, which 
acts as the prelude of planetesimal formation. The 
particle clumping generally requires the presence 
of relatively large particles with r s > 10~ 2 . It also 
sensitively depends on solid abundance, in favor of 
super-solar mctallicity. 

4. The dense particle clumps are mostly made of 
the largest particles with size range spanning less 
than one order of magnitude. These particles are 
trapped in the clumps for several orbital times be- 
fore leaving the clumps, providing a way for large 
particles to survive radial drift. 

5. The mean radial drift velocity for each particle 
species agrees well with a multi-species NSH equi- 
librium solution (see Appendix A). Strongly cou- 
pled particles drift outward, and the radial drift 
velocity for particles with larger r s is strongly re- 
duced relative to the conventional single-species 
NSH value, especially at large Z. This can increase 
the lifetime of the largest particles by a factor of a 
few. 

6. Turbulence generated by the SI leads to radial dif- 
fusion of particles, but the diffusion is slow and its 
effect is negligible compared with radial drift after 
about 100 orbital periods (for the largest particles). 
Particle clumping effectively enhances radial diffu- 
sion by retaining a fraction of large particles in the 
clumps. 



Fig. 11. — Summary of logical relations between various factors 
relevant to planetesimal formation. A solid arrow from A to B indi- 
cates that A promotes B. Dashed double arrow indicates speculated 
connection between the two processes. Z lar s e represents metallic- 
ity for large particles (millimeter or larger), Z sn represents the 
abundance for micron-sized small grains. Radial pile-up denotes 
the pile-up of particles due to the positive radial gradient of the 
radial drift velocity. We refer to "Grain growth" as collisional co- 
agulation towards the largest solid size. "Dead zone" refers to the 
radial extent of the dead zone. See text for details. 

7. Mutual collision velocity between t s > 10~ 2 par- 
ticles is dominated by the difference in their ra- 
dial drift velocities, and agrees well with calcula- 
tions using the multi-species NSH equilibrium. The 
collision velocity is strongly reduced towards large 
disk metallicity relative to predictions from single- 
species NSH solution. Collision velocity induced by 
SI turbulence is only secondary. 

7.2. Implications for Planetesimal Formation 

In this subsection we combine the results summarized 
in the previous subsection and discuss various implica- 
tions for planetesimal formation. In particular, we em- 
phasize that the importance of the local enrichment of 
solid materials in PPDs on planetesimal formation (two 
feedback loops, sec §7.2.21 and §7.2.3j) . Our logical chain 
is summarized in Figure 111! and we elaborate various 
aspects of this diagram in the following. 

7.2.1. Conditions for Strong Particle Clumping 

Our simulations show a dichotomy in the parameter 
space in which strong particle clumping occurs. Strong 
clumping requires the presence of relatively large parti- 
cles with r s > 10 _2 , and the abundance of these parti- 
cles to be super-solar. These two requirements are rep- 
resented by the two arrows connecting Z lal ' 6C and "grain 
growth" to "planetesimals" in Figure [TTJ These two re- 
quirements can compensate each other: to form plan- 
etesimals, less grain growth is required if the solid abun- 
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dance is large enough. In the case that all solids are 
strongly coupled to the gas, enhancing the disk metallic- 
ity by a factor of 5-30 may be able to trigger GI followed 
by planetesimal fo rmation (|Sekival[J998t lYoudin fc Shul 
|200|[Chian3[20Ql). 

The particle size that corresponds to r a = 10~ 2 de- 
pends on the disk model and distance from the central 
star. In the MMSN model, this stopping time corre- 
sponds to 2cm particles at 1AU, and 2mm or smaller 
particles at 5AU or beyond, according to equation ©. 
These particle sizes are about the maximum particle 
size obtained in th e recent dust coagulation calculations 
(jZsom et al.ll201(J ). For more massive disks, the corre- 
sponding particle size will be enlarged by a factor of f g 
in the Epstein regime. Since grain growth becomes more 
difficult when the particle size exceeds a millimeter, the 
SI scenario of planetesimal formation prefers less massive 
disks, or the outer part of the disk. 

7.2.2. Enrichment of Local Solid Abundance 

Enrichment of the local abundance of solids is possible 
by several effects, already briefly mentio ned in jfl] Here 
we foc us on the mechanism proposed by lYoudin fe Shul 
(|2002h . For a particle at a fixed size (small, with t s < 1 at 
all disk radii considered), the (single-species) NSH radial 
drift velocity in PPDs decreases as the particle drifts in- 
ward. Therefore, radial drift causes the particles to pile- 
up towards the inner regions, leading to enhancement of 
the local abundance of solids. Depending on the time 
scale, the largest enhancement factor can reach 3 — 10 in 
10 5 - 10 6 years (|Youdin fc Shull200J lYoudin fc Chiang 
12004 ) . starting in the outer disk and moving inwards. 
This process corresponds to the arrow pointing from "ra- 
dial pile-up" to Z Iargc in Figure HU 

In this paper, we have shown that the radial drift veloc- 
ity is reduced when local abundance of solids increases. 
This effect provides positive feedback to the enrichment 
process: particles pile up because the radial drift velocity 
is smaller at smaller disk radii. The enrichment of solid 
material at small disk radii further reduces the radial 
drift velocity, leading to stronger pile-ups. This effect 
corresponds to the arrow pointing from Z lar s c ^ " rac Ji a i 
pile-up" in Figure [TTJ We see that the enrichment of 
the abundance of solids and the particle pile-up form a 
feedback loop that enhance each other. Therefore, we 
expect even stronger solid enrichment i n the inner re- 
gion of PPDs than previous c alculations (|Youdin fc Shul 
l200l lYoudin fc Chiandl200l sufficient for SI and/or GI 
to form planctcsimals within PPD dead zone. 

7.2.3. Implications for Grain Growth 

The radial drift velocity adopte d in grain coagulation 
models (e.g.. IBrauer et alJ l2008at iBirnstiel et"atll2010t 
IZsom et al.ll2010D is generally taken from a single-species 
NSH equilibrium. However, we have shown that as par- 
ticles settle to disk midplane, the radial drift velocity is 
reduced due to multi-species effects, and smallest parti- 
cles can even drift outward rather than inward. More 
sophisticated modeling of grain growth is needed to in- 
corporate the multi-species effects. 

One consequence of the multi-species NSH equilibrium 
is that the enhancement of local abundance of solids 



strongly reduces the radial drift velocity, hence the par- 
ticle mutual collision velocity. Because the radial drift 
velocity dominates the collision velocities for relatively 
large particles, we expect particles to grow larger in re- 
gions with large abundance of solids, due to the reduced 
collision velocity there. This effect is illustrated as an ar- 
row pointing from Z lar s c to "grain growth" in Figure [TTJ 
In turn, grain growth form large particles from smaller 
ones, increasing the population of large particles, hence 
Z large , leading to the arrow pointing in the opposite di- 
rection. Again, these two effects form a feedback loop, 
promoting grain growth and enrichment of solid material. 

7.2.4. Influence from Magnetic Activity 

All our simulations ignore external sources of turbu- 
lence, particularly the MRI, by working in the dead 
zone of PPDs. Even very weak external turbulence with 
a ~ 10~ 4 may stir up the solids and maintain them at a 
height that may be insuffi cient for SI to form den se parti- 
cle clumps (see Table [2]). iJohansen et al.l (|2007|) showed 
that planetesimal formation is facilitated by MRI due 
to the long-lived overdensity regions that effectively trap 
the particles. However, they chose very large particle size 
with t s = 0.25 — 1. When smaller particles are used, they 
are diffused to a much larger height. Although SI may 
still be pre sent with MRI turb ulence and form elongated 
structures (|Balsara et al.l2009l ) , particle overdensities are 
small. Moreover, the mutual par ticle collision velocity i s 
much higher in MRI turbulence (jCarballido et al . 1 120081) . 
which inhibits particle growth. Therefore, we expect the 
dead zone to be the more favored region for planetesimal 
formation, and planetesimal formation should be easier 
in PPDs with a larger (radial) extent of the dead zone. 
This is shown as the arrow pointing from "dead zone" 
towards "planetesimal formation" in Figure [TTJ 

At the same time, the vertical and radial extent of the 
dead zone in PPDs strongly depe nd on the abundance of 
micron-sized and smal ler grains (jllgner fc Nelsonl 120061 : 
iBai fc Goodma n 2009), which is reflected in the arrow 
connecting Z sma11 and "dead zone". Whether there is 
any connection between Z largG and Z smldl is uncertain. 
Both Z lar s° and Z smaI1 should be correlated with the 
overall disk metallicity. Moreover, with more large parti- 
cles (larger 2" largc , collisional fragmentation may lead to 
more small grains (larger Z smaU ). Despite many uncer- 
tainties, we draw a dashed double-arrow between Z large 
and Z sma11 as our speculated connection between the two 
particle populations. If such connection exist, it repre- 
sents the third way for (local) solid abundance enrich- 
ment to promote planetesimal formation. 

7.3. Limitations and Outlook 

Our simulations take a local patch from a sim- 
ple global disk model in which all physical quantities 
follow a power-law dependence on disk radii. The 
global structure of PPDs may be more complicated. 
In particular, the presence of a dead zone in PPDs 
changes the steady-state disk surface density profile, and 
may lead to local pressure maxima at the snow line 
(jKretke fc Linl 120071 : IBrauer et all l2008bD. and the in- 
ner edge of the dead zone ( Dzvurkevich et all |2010[) . 
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These pressure bumps are able to trap particles very effi- 
ciently. M oreover, solids can al s o be t r apped in long-lived 
vorti c es (iBarge fe Sommerial 119951: iKlahr fe Henninsl 
19971: IKlahr fe Bodenheimerl 120031: iJohansen et al.ll2004l: 
Lvra et all l2009f ). although their existence is in de- 
bate. Finally, the structure of PPDs ma y be very 
non-s teady and undergo periods of outburst (Zh u et all 
I2009T) . Global models take into account large-scale varia- 
tions of disk structure and can follow the disk evolution. 
Currently, the numerical resolution required for resolving 
the SI is prohibitively high for running 3D global simu- 
lations. Nevertheless, one can perform local simulations 
in different disk environments , and piece them t ogether 
to form a global picture, as in iTillev et al.l (|2009f ). 

Our simulations focus on the dynamics in the vicin- 
ity of the disk midplane with very limited radial 
and vertical box sizes and no magnetic field. This 
is mainly constrained by the fine grid resolution re- 
quired for this study. In reality, MM turbulence in 
the active lay ers may excite vertical oscillations in the 
disk midplane ([Fleming fc Stonell2C)03l: lOishi fc Mac Low! 
|2009|) . Moreover, turbulent diffusion of ions may produce 
magnet ic activity at the dis k midplane, making it "un- 
dead" (jTurner fc Sanoll2008t ). Finally, turbulent mixing 
of particles may be come important wh en the active layer 
is relatively thick ([Turner et al.l l2010h . Including these 



effects into numerical simulations of planctesimal forma- 
tion requires enlarging our box size in all dimensions by 
a factor of at least 10. Moreover, 3D rather than 2D sim- 
ulations are necessary to maintain sustained MRI turbu- 
lence. Such simulations are computationally expensive, 
however recently the static mesh refinement (SMR) al- 
gorithm in the Athena MHD code has been parallelized. 
The cost of hybrid (particles and gas) simulations of lay- 
ered disks can be substantially reduced by using fine res- 
olution at the disk midplane (to capture the SI) and in 
the active layers (to capture the MRI) , while using coarse 
resolution everywhere else. With SMR, it becomes fea- 
sible to study the effect of (non-ideal) MRI turbulence 
in the active layer on particle dynamics and planetesi- 
mal formation in the disk midplane. This is planned as 
future work. 
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APPENDIX 

A. MULTI-SPECIES NSH EQUILIBRIUM 

In this appendix we generalize the NSH equilibrium solution to include multi-species of particles. We start from the 
force balance for both gas and particle components. Their velocities are denoted by u for gas and v k for particle type 
fc. Tfc and e/j denote dimensionless particle stopping time and solid to gas mass ratio for particle type fc. In writing 
down the equations, we subtract both gas and particle velocities by linear Keplerian shear — (3/2)Qxy, and denote the 
remaining velocities by u' and v' k respectively. Neglecting vertical gravity, the hydrostatic equilibrium equations read 



2v 'k v x ~ W kx y K - u') = , 

2 Tfc 



1 



2u' y x - -u' x y 



y] —( v 'k - «') = -2vvkx 



(Ala) 
(Alb) 



Multiplying equation (|Ala[) by e k for each fc and adding them to equation (|Alb[) . we find the expression of the gas 
velocity in terms of particle velocities 

(A2) 



k 



To obtain the particle velocities, we define velocity vectors T x = (v[ 



xi u 2xi 



,) T , and T y = (v[ 



u 2y i ■ ■ ■ i v n 



I \T 



Further we define a diagonal matrix A = diag-jYi, T2, . . . , r„}, and a rank 1 matrix T = (e, e, . . . , e) , where e 
(ei, €2, ■ ■ ■ , e n ) T - With these notations, the equations governing v' k can be written as 



fl + T -2A 
^ A/2 1 + r 



T x 

T„ 



= -TjVK 



(A3) 



The solution of this equation can not be expressed simply; but taking advantage of the block structure of the matrices, 
one can find the solution in the form 



T, 
T„ 



-rjVK 



A IB 
■B/2 A 



where 



s = {[A- 1 (i + r)] 2 + i}- 1 A- 



A = A- 1 (1 + T)B . 



(A4) 



(A5) 



One can easily verify that equation (|A4[) indeed generalizes the single species NSH solution. The multi-species 
solution obtained here is useful for setting initial conditions of the simulation as well as analysis of simulation data. 
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Fig. 12. — The cumulative probability distribution function (CPDF) of two-body particle collision velocities from runs R21Z3 (top), 
R30Z3 (middle) and R10Z3 (bottom), illustrating the effect of particle clumping on collision velocities. Results from both 2D (solid) and 
3D (dash-dotted) simulations are shown. In each panel, we plot the CPDFs for collisions between various particle species (with different 
colors, sec figure legends) and one specific particle species with fixed t s (indicated in each panel). Dashed lines mark the median value, 
while dotted lines mark the la level of fluctuations. Note that 3D runs in R21 and R30 do not show particle clumping, while the others do. 



B. COLLISION VELOCITY DISTRIBUTION 

In this appendix we discuss the distribution of particle collision velocities. In Figure IT21 we show the CPDFs from 
runs R21Z3, R30Z3 and R10Z3 (for both 2D and 3D). Each panel plots the collision velocity CPDFs between several 
(or all) particle species and a given particle species. The median of the collision velocity is represented by the dashed 
lines at P = 0.5, with la range obtained by cutting the plots at P = 0.68 and P = 0.32, shown in dotted lines. 

Two 3D runs R21Z3-3D and R30Z3-3D do not show strong particle clumping, and their CPDFs are representative 
of runs without particle clumping. The CPDF curves for different pairs of particle types are very similar between 
each other. The curves are generally anti-symmetric with respect to the median value P = 0.5, and the corresponding 
velocity distribution is close to log-normal distribution with super-exponential cutoff at large velocity. Moreover, even 
in runs that show strong particle clumping, the CPDFs for collisions between relatively small particles with t s < 10 -2 
also approaches the log-normal form, as one can see from run R30Z3-2D. The la limit of the collision velocity is 
generally less than half of the median collision velocity. 

The CPDFs from simulations that show strong clumping deviates substantially from log-normal. For runs R21Z3 
and R30Z3, the collision velocity is clearly reduced in 2D relative to 3D (2D runs show clumping while 3D runs do 
not). In run R30Z3-2D, we see that the reduction is most significant for collisions between two relatively large particles 
with t s > 0.1, which are the ones that make up most of the clumps (see Figure [6]). The collision velocity between large 
(t s > 0.1) and small (r s < 0.1) particles is also reduced, although to a lesser extent. For R21Z3-2D and R10Z3 (both 
2D and 3D) runs, all particles actively participate in the SI and clumping, and their mutual collision velocities are all 
reduced. The reduction of collision velocity appears to be stronger in R30Z3 and R10Z3 runs, where T max = 1. Finally, 
we see that the reduction of collision velocity is most significant for intermediate high- velocity collisions (V co i — 1 — 5m 
s _1 ), although there is still a high-velocity tail present. It is likely that the high-velocity tail is caused by collisions 
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Fig. 13. — The mean particle collision velocity as a function of particle density for the three runs that show particle clumping (R21Z3, 
R30Z3 and R10Z3, from top to bottom). The 2D results are shown in solid lines, while 3D results are shown as dash-dotted lines (bottom 
pan els only). In each panel, we show the median collision velocity between various particle species (using the same color scheme in Figure 
1 121 ) and a specific particle species with a t s as indicated in each panel. 

outside the particle clump. 

To further address the influence of clumping on particle collision velocities, we evaluate the mean collision velocity 
as a function of ambient particle density p p . For the three sets of runs that show particle clumping R21Z3 (2D), R30Z3 
(2D) and R103 (2D and 3D), we show the results in Figure [L3l For both R30Z3 and R10Z3 runs, we find a clear trend 
that the collision velocity is strongly reduced towards higher p p . More specifically, the reduction is most prominent 
for collisions between large and small particles. The collision velocity between particles of similar sizes appears to be 
insensitive to p p , and is maintained at a relatively low value. Interestingly, results from run R21Z3-2D show a different 
behavior. The collision velocity decreases with p p until p v j p 9i b ~ 300. Beyond this density, we observe an increase 
of collision velocity towards larger p p . We have performed an additional R41 run in 2D with Z=0.05, which shows 
particle clumping, and the collision velocity shows similar properties as in R21Z3-2D. This is very likely to be due to 
the different properties of the SI for p articles with r s = 0.1 from particles with r 8 = 1. As shown in non-stratified 
simulations (jJohansen fc Youdinl 120071 ) . clumps are much more dynamic in the former case, which may lead to larger 
collision velocities. 
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